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Technical  Objectives 

This  grant  has  as  its  goal  the  development  of  a  theory  of  vortex  breakdown  capable  of 
identifying  the  physical  mechanisms  responsible  for  the  phenomenon,  and  upon  which 
may  be  built  practical  methods  of  prediction  and  control  of  the  phenomenon  and  its 
contribution  to  the  forces  and  moments  on  aircraft  and  missiles.  The  two-year  grant 
concluding  May  31,  1989  initiated  a  systematic  approach  towards  this  ultimate  goal  in  a 
multi-part  program  with  more  limited  objectives. 

•  The  first  pan.  Task  A,  sought  the  construction  of  fully  nonlinear 
axisymmetric  solitary  waves  in  vortices.  This  objective  was  identified  in  the 
expectation  that  the  strong  deceleration  characteristic  of  vonex  breakdown 
can  be  associated  with  large  amplitude  axisymmetric  waves.  This 
association  provides  a  clear  mechanism  for  the  dominant  features  of  vonex 
breakdown,  and,  since  such  waves  were  expected  to  be  amenable  to  a 
relatively  simple  and  compact  mathematical  representation,  this  provided  the 
prospect  of  a  kernel  for  rapid  algorithms  suitable  for  design  or  control 
purposes. 

•  The  second  part.  Task  B,  sought  ways  to  characterize  the  stability  of  the 
fully  nonlinear  solitary  waves  found  in  Task  A  to  three-dimensional 
perturbations.  It  has  been  suggested  (and  laboratory  experiments  support 
the  suggestion)  that  these  large  amplitude  axisymmtric  waves  become 
unstable  when  they  grow  to  some  critical  amplitude,  and  that  this  instability 
provides  the  mechanism  to  fix  the  equilibrium  position  and  strength  of  a 
vortex  breakdown  structure. 

•  The  third  part.  Task  C,  constituted  a  study  of  the  linear  stability  of  a 
mathematical  model  (due  to  Hall  and  Stewartson)  of  the  leading  edge 
vortex.  The  model  leading  edge  vortex  is  explicitly  known,  and  the 
objective  of  this  part  is  to  develop  a  simplified  way  to  describe  important 
instability  mechanisms.  The  hope  is  that,  once  Task  B  is  completed,  simple 
descriptions,  related  to  Task  C,  can  be  fined  to  the  instabilities  of  the  large 
solitary  waves  that  we  expect  to  find.  Simplification  is  expected  to  be 
important  in  designing  rapid  algorithms  for  design  or  control  purposes. 
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•  Because  the  general  level  of  scientific  understanding  of  nonlinear  transitions 
in  vortex  flows  having  a  component  of  velocity  parallel  to  the  vortex  axis, 
as  is  the  case  in  aerodynamic  vortices,  is  virtually  nonexistent,  we 
undertook  an  additional  study  (Task  D)  on  transitions  in  rotating  Poiseuille 
flow.  This  flow,  an  exact  solution  of  the  Navier-Stokes  equations, 
provides  a  convenient  example  to  discover  some  of  the  important  nonlinear 
dynamical  features  of  vortices. 

Other  aspects  governing  the  physics  of  vortex  breakdown  are  of  key  importance,  but  were 
outside  the  scope  of  this  two-year  grant. 


Status  of  the  Research 

Task  A  has  been  completed.  Preliminary  reports  were  given  in  [1]  and  [2],  and  a  paper 
reporting  the  full  details  has  been  submined  for  publication  [3].  The  manuscript  for  [3]  is 
enclosed  and  is  made  part  of  this  report  by  reference.  The  principal  findings  are  that  fully 
nonlinear  solitary  waves  do  exist,  and  that  their  properties  are  approximately  captured  by  a 
remarkably  simple  ansatz.  The  paper  also  discusses  large  amplitude  periodic  wave  trains, 
nonlinear  deformations  of  vortex  flows,  and  the  interrelationships  between  these  various 
classes  of  vortex  flows.  A  more  limited  mathematical  method  also  has  been  applied  to  this 
problem.  A  draft  paper  on  this  work  exists  ([4]):  completion  of  this  paper  was 
consciously  delayed  partly  because  numerical  work  (since  done)  was  thought  desirable  to 
make  the  findings  of  [4]  more  concrete,  but  more  importantly  because  the  work  [3]  was 
given  a  higher  priority.  We  plan  to  finish  and  submit  [4]  for  publication. 

Task  B  is  in  progress.  The  mathematical  problem  to  be  solved  here  is  a  linear  eigenvalue 
problem  for  a  set  of  partial  differential  equations  in  three  space  dimensions.  The  problem  is 
not  separable,  so  the  standard  normal  mode  analysis  technique  of  hydrodynamic  stability 
theory  does  not  apply.  Because  the  numerical  problem  to  be  solved  is  extremely  large,  an 
attack  on  the  problem  capable  of  resolving  the  physically  important  scales  required  the 
development  of  new  numerical  techniques  to  solve  large  algebraic  eigenvalue  problems. 

The  lack  of  efficient  methods  to  accomplish  this  has  been  a  stumbling  block  in  the  past. 
Under  Task  B,  such  methods  have  been  developed,  and  a  report  on  the  algorithms 
developed  is  given  in  [5],  which  is  also  enclosed  and  incorporated  in  this  report  by 
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reference.  Results  using  these  numerical  techniques  on  the  stability  problem  were  not 
available  by  the  termination  of  the  reporting  period. 

Task  C  is  partially  completed.  Our  work  on  the  Hall-Stewartson  model  of  the  leading  edge 
vortex  has  addressed  the  large  Reynolds  number  limit,  for  which  we  have  detailed  the 
existence  of  unstable  modes.  This  work  was  first  reported  in  [6].  A  draft  of  a  paper  [7] 
with  a  full  description  has  been  prepared.  We  have  not  yet  found  the  critical  Reynolds 
number  for  the  onset  of  instability,  since  this  requires  a  more  extensive  numerical  effort 

Task  D  is  partially  completed.  One  paper  describing  mathematical  aspects  of  the  rotating 
Poiseuille  (pipe)  flow  exists  in  draft  form  ([8])  and  will  soon  be  submitted  for  publication. 
Other  work,  detailing  the  nature  of  transitions  and  mode  interactions  has  been  done  and  has 
been  briefly  reported  ([9]).  The  latter  work  traces  a  series  of  transitions  in  rotating  pipe 
flows  to  chaotic  states,  and  a  paper  will  be  written  to  report  that  investigation.  We  also 
have  started,  but  not  yet  completed,  work  on  three  wave  resonant  interactions  in  rotating 
pipe  flows. 
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Abstract 

Large-amplitude  axisymmetric  waves  on  columnar  vortices,  thought  to  be  related  to  flow 
structures  observed  in  vortex  breakdown,  are  found  as  static  bifurcations  of  the  Bragg- 
Hawthome  equation,  equivalent  to  the  steady,  axisymmetric,  Euler  equations.  Non-tnviaj 
solution  branches  bifurcate  as  the  swirl  ratio  (the  ratio  of  azimuthal  to  axial  velocity.) 
changes,  and  are  followed  into  strongly  nonlinear  regimes  using  a  numerical  continuation 
method.  Four  types  of  solutions  are  found:  multiple  columnar  solutions,  corresponding  to 
Benjamin's  “conjugate  flows",  with  subcritical-supercritical  pairing  of  wave  characteristics: 
solitary  waves,  extending  previously  known  weakly  nonlinear  solutions  to  amplitudes  large 
enough  to  produce  flow  reversals  similar  to  the  breakdown  transition:  periodic  wavetrains: 
and  solitary  waves  superimposed  on  the  conjugate  flow  that  emerge  from  the  periodic 
wavetrain  as  the  wavelength  or  amplitude  becomes  sufficiently  large.  Weakly  nonlinear 
soiiton  solutions  are  found  to  be  accurate  even  when  the  perturbations  they  cause  are  fairly 
stronti. 


1.  Introduction 

This  paper  is  concerned  with  axially  symmetric  standing  wavetrains  ana  solitary 
waves,  without  restriction  to  mtinitesimai  or  weakly  nonlinear  amplitudes,  in  inviscia 
incompressible  vortex  tlows  .  Although  the  paper  may  be  regarded  strictly  as  a 
contribution  to  understanaing  of  waves  which  may  propagate  on  vortex  cores,  our 
motivation  is  the  exploration  of  a  conceptual  picture  of  vortex  breakdown  given  by 
Leibovich  (1983)  (and  in  a  more  widely  accessible  review  in  1984,  we  shall  designate 
either  of  these  references  as  L). 

In  aerodynamic  contexts,  the  global  flowtield  causes  impressed  forcing  on 
concentrated  vortices  embedded  in  it.  Generally,  the  spatial  scales  of  the  forcing  are  large 
compared  to  scales  associated  with  the  vortex  core.  A  conceptual  model  of  vortex 
breakdown  is  promulgated  in  L,  guided  in  pan  by  laboratory  experiments  and  in  pan  bv 
the  weakly  nonlinear  "trapped  wave"'  theory  of  Randall  and  Leibovich  (1973;  or  RL).  In 
the  scenario  outlined  in  L.  vortex  breakdown  is  a  process  that  involves  a  crucial  admixture 
of  a  strongly  nonlinear  axisvmmetnc  wave  propagating  in  a  vortex  "waveguide"  having 
axially  varying  cnaractensucs  tana  hence  an  axial  pressure  gradient),  and  a  smaller 
asymmetric  perturbation  arising  from  instability  of  the  big  wave.  The  pressure  gradient 
impressed  by  the  waveguide  does  work  on  the  axially  symmetric  wave,  causing  it  to  grow 
to  large  amplitude.  The  weakly  nonlinear  trapped  wave  theory  of  RL  indicates  that  w  ave 
growth  of  this  sort  leads  to  a  positional  instability  of  the  wave,  as  it  grows,  it  propagates 
faster,  and  an  equilibrium  position  cannot  be  established  (with  respect  to  a  reference  frame 
fixed  by  the  waveguide)  unless  there  is  a  mechanism  tor  extraction  of  energy  from  the 
axisvmmetnc  wave.  RL  invoked  viscosity  as  a  mechanism  to  dissipate  wave  energy.  On 
the  other  hand,  laboratory  experiments  (such  as  those  done  by  Sarpkava,  1972  and  by 
Faler  and  Leibovich.  1977,  1978)  show  clear  evidence  of  nonaxisymmetric  features  within 
a  nearly  axisvmmetnc  "bubble"  form  of  vortex  breakdown,  and  the  onset  of  asymmetry  is 
consistent  with  an  instability'  of  the  axially  symmetric  flow.  If  the  "bubble"  were  regarded 
as  a  manifestation  of  a  large  amplitude,  nearly  axisvmmetnc.  wave,  then  instability  to 
asymmetric  perturbations  offers  the  possibility  of  much  larger  energy  transfers  from  the 
wave  than  would  viscous  dissipation.  Furthermore,  the  features  of  vortex  breakdown 
appear  to  depend  little  on  viscosity,  at  least  at  higher  values  of  the  Reynolds  number. 
Consequently,  it  is  suggested  in  L  that  the  required  energy  extraction  from  the  strongly 
nonlinear  axially  symmetric  wave  arises  by  the  transfer  of  its  energy  to  azimuthaily 
asymmetric  modes  of  mouon  which  arise  by  instability,  tit  is  pointed  out  in  L  that  a  large 
amplitude  axially  symmetric  mode  -  large  amplitude  implying  a  variation  occuring  on  axial 
scales  comparable  to  the  vortex  core  radius  -  is  required  in  vortex  breakdown,  since  only 
this  component  of  a  Fourier  decomposition  in  the  azimuth  can  lead  to  the  deceleration  of 
fluid  on  the  vortex  axis  that  is  the  hallmark  of  vortex  breakdown.)  On  the  basis  of  their 
experimental  observations  in  flows  very  different  from  tnose  already  cited.  Maxworthy. 
Morv,  and  Hopfinger  1 1983)  also  suggest,  in  a  paper  published  in  the  same  1983  volume 
as  the  paper  by  Leibovich.  that  breakdown  is  associated  w  ith  loss  of  stability  of  large 
axially  symmetric  waves  to  nonaxisymmetric  perturbations. 


To  explore  tnis  suggested  process,  we  have  divided  it  into  eiements  which  at  a  iaier 
time  must  be  recombined.  The  first  element  is  the  large  amplitude  axiailv  symmetric  wave, 
and  the  the  aim  or  this  paper  is  to  oeveiop  a  better  understanding  or'  these  waves.  The  otner 
ingredients  of  the  hypothesis,  not  yet  considered,  are  loss  of  stability  to  asymmetric 
perturbations,  and  effects  engendered  by  axial  inhomogeneity  caused  by  the  global 
flowfieid. 

In  axially  symmetric  steady  rlow.  it  is  well  known  that  the  Euler  equauons  for 
inviscid  swirling  flow  may  be  reduced  to  a  single  elliptic  paruai  differential  equauon  for  the 
streamfunction.  This  equation  seems  to  have  been  first  discovered  by  Bragg  and 
Hawthorne  ( 1950),  and  in  this  paper  we  will  refer  to  it  by  their  names.  Following 
Leibovich,  1985,  we  parameterize  the  Bragg-Hawthome  equation  (BHE)  by  the  relative 
level  of  the  swirl,  and  construct  standing  waves,  either  infinite  wavetrains  or  solitary 
waves,  by  studying  the  branching  behaviors  which  arise  as  the  parameter  is  varied  (other, 
nonwavv  flows,  were  discussed  by  Leibovich.  1985,  from  this  starting  point,  and  a  few  of 
the  results  to  be  detailed  here  for  weakly  nonlinear  wavy  flows  were  announced  in 
Leibovich.  1987).  The  procedure  can  therefore  be  described  as  a  search  for  static 
bifurcations  of  the  BHE.  Another  natural  parameterization  that  might  have  been  chosen  as 
a  starting  point  for  a  bifurcation  analysis  is  the  wave  speed  of  waves  of  permanent  form. 
Here  one  adds  a  constant  parameter  to  the  primary  axial  velocity  profile,  given  in  a  specific 
frame  of  reference,  and  regards  the  swirl  as  fixed.  This  choice,  although  more  useful  for 
some  purposes  than  the  parameterization  using  swiri  level,  introduces  the  wave  speed 
parameter  into  the  problem  in  an  awkwardly  nonlinear  way.  and  we  have  not  made  use  of  it 
in  this  work. 

Branching,  when  it  occurs,  is  from  a  primary  columnar  vortex,  assumed  to  be 
given.  The  flows  which  bifurcate  from  this  vortex  are  required  to  have  the  same  volume 
flow  rate,  and  the  same  total  head  and  circulation  variation  with  streamfunction  as  does  the 
given  vortex.  With  these  constraints  and  the  assumption  that  the  flow  is  either  periodic  in 
the  axiai  direction  with  a  finite  waveiengtn  L  or  is  columnar  at  upstream  and  downstream 
infinity,  we  find  that  new  flow  branches  may  be  of  four  kinds.  One  ciass  (I.  discussed  in 
H)  of  bifurcating  flows  is  again  coiumnar.  so  there  are  no  axiai  variations:  a  second  ciass 
( II.  §6)  consists  of  solitary  waves  with  the  primary  flow  at  upstream  infinity;  periodic 
wavetrains  comprise  the  third  ciass  (III.  §5);  and  a  fourth  class  (IV.  §5)  consists  of  solitary 
waves  that  approach  a  columnar  flow  at  large  axial  distances  that  is  disnnct  from  tne 
primary  columnar  flow. 

The  coiumnar  branches,  of  which  there  is  an  infinite  number,  u  hen  taken  together 
with  the  primary  flow,  are  the  "conjugate"  flows  defined  and  discussed  in  Benjamin  s 
( 1962)  seminal  paper.  Two  of  these  are  especially  important  and  receive  emphasis  in 
Benjamin  s  work:  these  are  the  primary  flow,  and  what  we  shall  call  the  principal  conjugate 
flow,  which  is  the  coiumnar  flow'  crunching  at  the  principal,  or  lowest,  eigenvalue  of  the 
linearized  problem,  and  therefore  differing  from  the  primary  flow  by  an  azimuthal  voracity 
(hat  does  not  change  sign.  The  procedure  adopted  here  allows  one  to  determine  how  these 
flows  are  connected  together.  We  are  also  able  to  show  that,  at  a  columnar-columnar 
bifurcation  point  between  the  primary  flow  and  its  principal  conjugate,  there  is  a  transfer  of 


Beniamin  s  1 1962)  criticality  classification  of  the  tiows.  That  is.  it' one  rlow  is  supercnucai 
on  one  side  of  the  bifurcauon  point  (does  not  admit  upstream  propagating  waves  or 
infinitesimal  amplitude,  and  shorter  standing  waves),  then  tne  other  columnar  orancn  is 
mtxnncai  (admits  upstream  propaganng  waves),  and  that  these  properties  of  the  branches 
.ire  exchanged  as  the  bifurcation  point  is  passed.  (N.B.  We  use  some  of  the  ianguage  of 
bifurcation  theory,  but  our  use  of  the  terms  supercritical'  and  suDcnticai  is  according  to 
Benjamin  s  wave  classification,  and  is  not  reiated  to  the  direction  of  bifurcauon.) 

Fully  nonlinear  solutions  have  been  tietermineu  numerically.  This  of  course 
requires  that  specific  examples  of  primary  voraces  be  considered.  We  have  chosen  one 
family  of  examples,  uniform  axial  velocity  and  the  Burgers  vortex.  While  the  uniform  axial 
velocity  of  this  primary  vortex  differs  from  the  jet-like  or  wake-like  character  of  flows  in 
which  vortex  breakdowns  have  been  observed,  the  principal  columnar  tiows  w  hich 
bifurcate  from  our  example  are  either  jet-like  or  wake-like  depending  on  swiri  level,  and 
therefore  closely  resemble  flows  upstream  or  downstream,  respectively,  of  vortex 
breakdown.  In  addition  to  the  general  results  on  columnar  branches  already  mentioned, 
we  note  here  an  interesting  -  thougn  academic  -  feature  of  the  principal  columnar  brancn. 

We  have  found  that  this  branch  can  be  numerically  continued  to  very  small  swiri  levels. 
With  this  clue  as  a  point  of  departure,  we  show  (in  $4.4)  the  limit  of  vanishing  swirl  by 
asvmptouc  means,  showing  that  connnuation  to  zero  swiri  is  possible. 

Numerical  examples  of  fully  nonlinear  wavy  solutions  similar  to  some  of  ihose 
discussed  in  this  paper  have  been  presented  by  Hafez  et  al.  (1986)  (an  abbreviated  account 
is  given  by  Hafez  and  Salas.  1985).  Our  paper  gives  a  considerably  more  comprehensive 
picture  of  the  inviscid  picture  by  identifying  classes  of  solutions  other  than  periodic 
wavetrains.  and  by  uncovering  the  connections  between  them.  The  results  of  Hafez  et  al. 
w’ere  computed  for  a  different  (although  more  or  less  similar)  class  of  primary  vortices,  yet 
their  results  are  qualitatively  similar  to  ours.  This,  and  the  theoretical  argument  of  this 
paper,  leads  us  to  believe  that  the  response  to  variation  of  swirl,  and  the  main 
charactensucs  of  the  tiows  to  be  described  here,  are  not  the  consequence  of  a  special  choice 
of  profiles  tor  the  numerical  worK.  but  are  of  general  applicability. 

Solitary  waves  of  class  II  can  be  found  for  the  supercritical  values  of  parameter  by 
numerically  continuing  known  weakly  nonlinear  solitar  waves  (with  our  choice  of  primary 
vortex,  we  have  access  to  weakly  nonlinear  results  found  by  Leibovicn.  1970.  w  hich  we 
use  as  a  starring  point)  to  large  amplitude.  These  solitary  wave  flowfields  approach  the 
primary  columnar  flow  asymptotically  at  large  upstream  and  downstream  distances.  Their 
amplitudes  can  be  increased  to  values  sufficiently  large  to  cause  stagnation  points  to  appear, 
followed  by  encapsulated  regions  of  closed  streamlines.  The  utility  of  the  results  when 
tlow  reversals  are  present  requires  careful  consideration,  particularly  since  added 
nonuniqueness  enters  in  such  circumstances,  and  the  consistency  requirement  ot  the 
Pranatl-Batchelor  criterion  is  violated.  We  believe  that  the  results  wiil  prove  to  be  of  value 
at  these  large  amplitudes,  and  discuss  the  reasons  for  this  belief  in  the  final  section  ($7). 

An  interesting  and  potentially  useful  point  is  that  even  at  amplitudes  large  enough  to  cause 
'tagnanon  and  reversed  flow,  the  weakly  nonlinear  solitary  wave  solutions  remains  an 
good  approximation  to  the  results  of  the  fully  nonlinear  calculations. 


For  suocnncai  parameter  values,  penouic  wavetrains  tciass  ill  )  have  oeen  found  by 
numerical  connnuation  beginning  with  infinitesimal  waves.  As  the  amplitudes  or  the 
periodic  wavetrains  increase  at  fixed  waveiengtn  in  the  numerical  examples  (as  they  uo 
when  the  parameter  measuring  swiri  increases),  the  wave  troughs  become  highly  localized, 
and  the  wave  crests  oecome  very  oroau.  These  oroad  crests  have  vinuatiy  no  axial 
■.  ariation.  that  is.  they  are  neariv  columnar.  We  show  tnat  tnese  nearly  columnar  rlow 
closely  approximate,  at  large  wave  amplitude,  the  columnar  principal  conjugate  vortex. 

The  same  kind  of  behavior  is  found  to  hold  when  the  swiri  level  is  held  fixed,  but  the 
wavelength  of  a  wavetrain  is  increased.  In  either  case,  an  individual  wave  trough 
approaches  a  solitary  wave  with  the  principal  conjugate  columnar  vortex  being  the  flow  at 
large  '  distances  upstream  and  downstream.  Thus,  we  are  led  to  another  family  of  solitary 
waves  (class  £ VI  distinct  from  those  supported  by  the  primary  vortex,  but  clearly  connected 
to  it. 


2.  Problem  formulation 

The  Euler  equations  in  cylindrical  coordinates  ir.B.z)  for  steady,  incompressible, 
axiaiiy-symmetnc  rlow  may  be  reduced  to  a  single  elliptic  partial  differential  equation  for 
the  Stokes  streamr'unction.  vj,  related  to  the  radial  (u)  and  axial  (wi  by 

u  =  -  r'  .  w  =  r*  ‘  y/r  . 

The  equation,  which  connects  the  azimuthal  voracity  to  the  total  head  and  circulation, 
seems  to  have  been  first  found  by  Bragg  and  Hawthorne  ( 1950)  and  is  given  by 


D“\y  =  r-H'(V)  -  EFriin 


BHE) 


where 


D~\j  =  vi 


rr 


I  / r )  'i/f  i-  *  • 


is  — r  the  azimuthal  component  of  voracity.  A  number  of  other  authors  have  made  use  of 
this  equation  notably  Long  1 1953).  Franxel  s  1956).  Squire  ( 19561.  and  Benjamin  t  1962). 
See  Batchelor  ( 1967)  for  a  convenient  reference  for  a  derivation  of  BHE.  following  Bragg 
and  Hawtnome.  or  Yih  ( 1965  )  for  a  derivation  by  an  alternate  method.  The  replacement  ot 
trie  Euler  equations  ny  BHE  is  justifiable  at  ail  points  in  a  steady,  inviscid  and  axially 
symmetric  flow  with  the  possible  exception  of  meridional  stagnation  points  ti.e..  wnere  u  = 
w  =  0.  but  v.  the  azimutnal  velocity  component,  may  he  nonzero),  which  are  singular 
points  of  the  transrormaiion.  In  BHE.  Him/),  the  totai  head  or  Bernoulli  function 


IFlI/l  =  —  n  -*-  <  i/2>  v»v  . 

P 

ind  FCij/),  the  circulation  aoout  ine  symmetry  axis  lapan  irom  a  factor  of  2rti. 


Ft'vi/)  =  rvir.zi  i2) 

are  functions  of  vi/  aione.  That  Fisa  funcuon  of  mi  aione  is  a  consequence  of  conservation 
of  anguiar  momentum  in  an  mviscid  fluid. 

Solutions  to  ( BHE)  are  solutions  to  the  Euler  equations  under  the  stated 
restrictions,  and  different  mviscid  flow  problems  arise  from  the  specificauons  of  the  pair  of 

functions  H  and  F.  The  simplest  cases  are  those  for  which  the  radial  component  of  velocity 
vanishes  at  some  plane  z  =  z.  on  which  the  axial  (w>  and  azimuthal  t\  >  velocity 

components  are  specified  to  be  W(r),  Vfr)  with  vV  ■*.  0.  then  HYun  and  F  are  easily 
determined  (see  Beniamin.  1962).  In  this  case,  which  we  consider  here. 


1  dp  V— ( r ) 


at  this  plane.  Since 


vj/(r,z, )  =  f  rWfr)  dr  nvF(r) 

‘  0 


we  can  suppose  this  latter  relation  to  be  invened  to  give 


r  =  Rim/) 


Fhis  inversion  can  aiwavs  be  done  un  pnnciple)  if  W(r)  t  i).  At  z  =  /  we  can  now 
regard  W  as  a  function  of  \j.  since 

!  dw  1  dM/  (0) 

w  =  t  -rr  ir-z,  >  = - 

1  ar  ‘  R(vt/)dK 

With  V<n  presented.  Fng)  i>>  determined  to  be 


:  V  =  Fi vri  =  Ruin  \'i R( \u )- 


Hie  pressure  may  now  be  found  as  a  function  of  M /  from  1 3 1  by  intetrration. 


R{  Vi/)  i.N) 

-N 

i  i  .  i  V  ~  ( r )  . 

—  p  =  —  mO.z.)  *  |  — - — ur. 

P  P  (1 


and  thus  H(V)  may  be  identified.  If  the  raaial  velocity  is  not  given  to  be  zero  at  z=Zj.  but  is 

prescribed  as  some  nontrivial  function  of  r,  a  similar  construction  ot  H(\j/),  F(Vi/)  can  be 
carried  out.  We  call  the  flow  given  at  the  piane  z=Zj  the  specifying  flow  .  this  wiil  be 

taken  to  be  the  basic,  or  primary,  flow  and  the  stoning  point  of  our  investigations. 

Only  H'lVj  is  required  for  the  analysis.  H(i|/)  itself  is  not  needed.  At  a  piane  z=z5 

where  u  =  0,  v  =  V.  w  =  W,  and 


dH  _  PHdR  =  1  fdF  1  dW 
dvi/  r)r  dvt/  R“  ..iur  ^  ^ 


Now  we  suppose  that,  at  z  =  z,,  the  functional  form  of  the  swirl  Vfr)  is  nxea.  but 
the  level  is  adjustable,  so  that 

Ffy)  =  a  f(g/)  (10) 

where  f  is  a  fixed  function,  and  X  is  an  adjustable  constant.  Equation  ( BHE)  fcr  <jr  iSfay 

-* 

be  written 

D“vit  =  rAlwi  -  X"B(g/,r")  (11) 

where 

A (\\n=— — ^iR(u/)),  and  Blti/.r")  02) 

R(\j/)  dr 

=  . R2(\|/)-r~)  ff  '(virt/R^d/)  . 

With  the  appropriate  interpretation  of  X.  we  may  regard  equation  1 1 1 )  as 
dimensionless.  Thus,  if  we  scale  distances  by  a  characteristic  radius  b  (such  as  that  ot  a 
bounding  tube.  or.  alternatively,  the  location  of  the  maximum  swirl  speed),  the  specitving 
axial  velocity  with  a  characteristic  speed  Wo  (such  as  its  value  on  the  axis),  (he 
>treamtuncuon  with  b-Wo.  the  azimuthal  speed  with  a  typical  value  V.  it  such  as  the 
maximum  occunng  in  tne  flowi  then 


Vo 
Wo  ' 


M- 


ana  we  may  interpret  i  i  1 )  as  a  uimensioniess  equation. 

The  parameter  /.  appears  ontv  as  k-.  and  so  we  re  mace  it  nereaiter  wan 


and  as  a  consequence  of  its  definition,  only  admit  posidve  values  of  A. 

Suppose  the  boundaries  of  the  fluid  in  the  (r.zi  plane  have,  as  two  constituents,  the 
impermeable  cylinders  r  =  a  and  r  =  1  >  a  (here  we  have  chosen  the  outer  rube  raciius  as 
length  scale  for  our  problem).  The  columnar  specifying  tlow 

\'J  =  vFfr)  .  ( 14) 


is  a  solution  of  ill)  holding  for  ail  A  >  0.  We  now  wish  to  find  other  solutions,  periodic 
m  z  with  a  prescribed  wavelength  L.  Since  both  the  specifying  columnar  ilow  and  any 
other  wave-iike  solutions  that  may  exist  simultaneously  satisfy  the  same  mathemancal 
problem. and  one  possibility  for  this  to  occur  is  by  bifurcation  of  new  branches  of  solutions 
from  the  specifying  tlow,  and  the  multiple  solutions  so  obtained,  when  of  small  amplitude, 
rnav  be  identified  as  the  waves  propagating  on  the  specifying  tlow  previously  found  m  the 
literature  tcf.  Long.  1953,  Frankel.  1956,  Squire.  1956.  Benjamin.  1962).  Lei  -  'f'- 


O  =  \j  -  vP(r), 


H5) 


represent  the  perturbation  streamfunction.  If  there  are  other  solutions,  there  is  a  z-penodic 
nontnviai  solution  to  the  (elliptic)  partial  differential  equation 


where 


MO.  A)  =  D“0  +  iiiO.r.A)  =  ■) 


Ll(0.r.k)  =  A  P(O.r)  -  r-Q(O.r) 

P(O.r)  =  B(VF  h-  O.  r“ )  -  B(VF.  r  > 
Q(O.r)  s  A(4,  +  <P)-A(H') 


•  16b) 


atisfvine  the  bounaar\'  conditiorw 


O(a.z)  =  0(  1  ,z)  =  (>. 


(I)(r,z  •  tLi  =  <l)(r.z  *tL) 


The  nonlinear  proolem  admits  solutions  even  in  z.  and  we  focus  on  tnese.  In 
addition  to  admitting  solutions  with  this  symmetry  iz  —  -  zi.  soiuuons  are  aiso  aarruttea 
with  z  — r  z  •+•  h.  tor  any  n.  Thus,  smooth  z-penodie  solutions  may  be  constructed  by 
appropriately  piecing  togetner  (by  reflections  and  shifts)  soiuuons  satisfying  the  Neumann 
(boundary  conditions. 


do  do  1 

—  r.U>  =  — T.tD=». 


3.  Static  bifurcation  analysis 

The  Bragg-Hawthome  equation  describes  only  steady,  or  static',  solutions.  It  can 
therefore  be  used  to  describe  branches  of  the  families  of  steady  solutions  corresponding  to 
the  same  funcuonai  forms  (as  functions  of  the  streamfuncnon )  for  the  total  head  and 
circulation,  and  the  same  voiume  rate  of  flow.  The  bifurcation  and  connnuauon  of  such 
branches  is  discussed  in  this  section.  The  quesdon  of  the  stability  of  the  various  solution 
branches  of  the  BHE  is  a  dynamical  problem.  This  cannot  be  answered  in  the  context  of 
the  BHE  equation,  and  it  is  necessary  to  return  to  the  Euler  equations  in  which  solutions  to 
BHE  are  embedded.  This  is  addressed  (for  columnar  solutions  only)  in  §5. 

3.1.  Perturbation  expansion 

We  know  that  the  specifying  flow.  =  0,  is  a  solution  to  the  problem  1 16)  for  any 
value  ot  the  parameter  A.  This  can  fail  to  be  a  unique  solution  branch  for  a  given  A  only 
when  L(O.A),  the  operator  defined  by  the  linearization  of  N(O.A)  about  O  =  0.  is  not 
invertible.  This  occurs  oniy  when  the  parameter  A  coincides  with  an  eigenvalue.  \i  t  say ). 
for  the  linearized  nrooiem 

L(O.A)00  =  D-"0n  -r  — ->0.r,A)  0Q  ‘ 17a) 

">  (/P  ~T)C) 

=  D“0n  -r  |  A — ’O.r)  -  r~ — (0,rl|Qn  =  0. 

TO) 

O0(a.zi  =  o()( l.z)  =  0.  onva.z  -  tL)  =  on(a,z  +|L).  (17b) 

Tor  A  near  an  eigenvalue  u.  we  construct  a  solution  branching  from  me 
'pecitying  flow  in  a  perturbation  senes.  This  wiil  provide  a  local  approximation  for  the 
olution  branching  at  u.  which  we  wiil  continue  numerically  to  iarcer  values  of  (A  -m. 


<t>  =  !'(0(,  -r  t'O, 


.\(E)  =  U  +  EK(£)  =  U  +  t'(K*0  -  EK,  ~ 


ana  set 


1  <5kP  M 

Pkir)  =  -T— 0-r). 
K-  dd>k 


;  1 8a ) 
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qktd  = 


i 

k!  aok 


<0.r). 


(19b) 


Sj.(r;f!)  =  upk(r)  -  r~qk(r),  ,19c) 

where  u  refers  to  any  eigenvalue  of  the  linearized  proDlem  1 171  ana  o<)  the  corresponding 
eigenfunction,  e  is  a  small  ordering  parameter  measuring  the  amplitude  of  the  bifurcating 
solutions  and  the  difference  between  A  and  its  value  p.  at  the  bifurcation  point,  and  the  dots 
stand  for  higher  order  terms  in  e.  When  these  relationships  are  substituted  into  (16a). 
Taylor  senes  expansions  in  powers  of  e  earned  out  and  each  coefficient  in  the  senes  is  set 
to  zero,  the  first  three  coefficients  are 


UO,p.)({>0  =  D200+  p,(r)  -  r~q,(r)|0o 

LTO.pjo,  =  -  (tc()p,(r)0o  ^  s2(r:)Li)Oo“ } 

UO.pnOi  =  -  !  K.pjtnOo  +  xypprio,  +  K*„p2tr:u)0(r  - 
2s2t  rigt  )0n0,  +  s -strip. kv  ! 


and  all  the  on  satisfy  the  same  boundary  conditions  1 17b). 
Note  that 

1  df“ 


p,(r)  s  - - — 

W“(n 


and 


M,(f) 


1  d  .  1  dW 
rW  dr  r  ar 


(20a) 
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The  numerator  or  p,tr),  when  muinpiied  by  A.  is  Rayleigh  s  discriminant,  ana 

therefore  sensible  orooiems.  in  the  context  of  considerations  such  as  in  this  paper,  wiii 
.have  p|(ri  >  0.  and  we  assume  this  to  be  the  case.  (Otherxvi.se.  the  primary  now  is 

unstable.)  Only  positive  values  of  u  can  correspona  to  brancn  points,  since  A  >  0  by 
definition.  If  p,irj  >  0  ana  q,(n  >  0.  then  the  smallest  (or  principal")  eigenvalue  is 
guaranteed  positive  (Leibovicn,  1985).  and  is  therefore  a  possible  branch  point.  Positive 
values  of  q  i  (.  r ;  are  not  necessary  for  this  to  be  so.  In  the  development  below,  we  tacitly 
assume  the  smallest  eigenvalue  is  positive.  If  this  is  not  so.  then  the  mathematical  changes 
needed  are  obvious  -  one  deals  oniv  with  the  positive  eigenvalues,  of  which  there  is  an 
infinite  number  -  but  the  branching  solutions  are  likely  to  be  unstable  and  therefore 
physically  uninteresung. 

The  problems  we  have  posed  here  depend  on  two  parameters.  A  and  L.  once  the 
specifying  flow  is  selected. 

The  pnncipai  eigemuncnon  has  no  zeros  in  the  intenor  or  D  fCourant  and  Hilben. 
1953),  and  without  loss  of  generality,  we  therefore  may  take  it  to  De  nonnegative.  Even 
eigenfunctions  are  all  of  the  form 

Oo  =  /.pmtnc°s(2nmz/L).  (21) 

where  the  integer  index  p  is  the  number  of  internal  zeros  of  Xpmir),  which  satisfies  the 
problem 


d 

dr '  r  ar 


2mn  2, 

mPpn-r-qpn  i  ZPm 


=  0 


(22) 


Xpm^*  /pin'1*  **• 


Here  we  nave  iabeled  the  eigenvalues  according  to  the  indices  p  and  m  corresponding 
to  the  associated  eigenfunction.  The  pnncipai  eigenfuncuon  corresponds  to  m  =()  and  to  the 
index  p,  which  we  can  take  to  be  p  =  0.  such  that  the  funcnon  £,K)(r)  has  no  zeros  internal 
to  D.  The  pnncipai  eigenfunction  belongs  to  the  eigenvalue  a,K)and  is  a  function  of  r  alone. 
A  solution  which  branches  from  the  pnncipai  eigenvalue  therefore  corresponds  to  a  new' 
columnar  flow  i  which  we  will  call  the  "principal  conjugate  branch"  since  it  is  a  conjugate 
(low  as  definea  by  Beniamin.  1962).  and  an  infinite  number  of  other  columnar  flows  i  also 
contugates  in  the  sense  of  Beniamin )  branch  from  iarger  eigenvalues  corresponditm  to  in  =  <• 

..nd  the  eigenfunctions  /putn.  for  p  =  i  .2.3 .  Solutions  periodic  in  z  (standing  waves i 

branch  from  eigenvalues  corresponding  to  eigenfunctions  w  ith  m  ^  0.  Modes  for  ail 
values  of  m  are  characterized  by  the  number  of  zeros  of  their  eigenfunctions  w  ith  7;im;n 
having  p  internal  zeros. 


;o- 


The  eigenvalue  problem  1 1 7 »  is  in  stunaaru  Sturm-Liouviile  torm  iCourant  a; 
Hilbert.  1953).  ana  some  ot'its  features  isuch  as  Dounas  on  me  smallest  eigenvalue;  are 
discussed  bv  Leibovicn  (1985).  There  is  one  comment  which  is  worth  maxing  at  this  point 
about  this  eigensvstem.  in  addition  to  the  observations  we  have  already  made. 

Eigenvalues  corresponding  to  z-depenaent  eigenfunctions  (constituting  w  avy  modes, 
with  eigenvalues  exceeding  fXxd.  decrease  as  L  increases  (Courant  <&  Hilbert.  1953 1.  ana 
as  L  -4  °o.  is  an  accumulation  point  tor  the  wavy  eigenvalues  m  r-  0.  For  tne 
->ame  reason,  waves  corresponding  to  higher  radial  modes  have  accumulation  points,  with 
dpm  M-po  us  L  — > 

3.2.  Branching  behavior 

If  do  is  an  eigensolution.  then  the  solution  to  the  adjoint  eigenvalue  problem 
with  an  unweighted  inner  product  is  Or/r,  or  alternatively,  the  problem  is  self-adjoint 
under  the  inner  product 

( F.G  i  =<  FG> 

where 

(•)  =  lr'*(»)drdz. 

D 

and  D  is  the  spatial  domain  in  which  our  problem  is  set. 

For  the  problem  for  Oj  to  have  a  solution,  the  solvability  condition 

TitniaiOo')  (33) 

Kq  —  —  "  i 

■  P|(riOn“) 


must  be  satisfied.  If  u  =  u,K)  is  the  smallest  eigenvalue,  then  o„  is  the  pnmary 

■v  ..... 

eigenfuncnon.  which  may  be  taken  to  be  positive.  L'nless  ipiop^m  -  r-u-qr;  is  aistnouteu 
in  a  special  way,  then,  x,  is  nonzero  and  the  bifurcation  at  (Xto  occurs  with  finite  siope 
i  i.e.,  it  is  transcruicai  or  dA/de  *  0  at  e  =  0). 

The  eigenfunction  corresponaing  to  j ,  the  lowest  eigenvalue  for  m  =  1.  is 

°o  =  Xn)tr)Cos(27tz/L).  (24) 

According  to  (23).  X,  =  0  for  solutions  branching  from  u,M,  once  the  defining 
integrals  extend  over  one  period  in  z  and  the  numerator  therefore  vanishes.  Hie  ditterenuai 
equation  determining  0(.  from  (20bt.  is  now 


i  1- 


25) 


L(0.(i)Oi  =  -  'M-oiP-yiri  -  r-q„in|  Of“ 


=  -  4i^„p,(n  -  r:q,(rj]| 
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with  solutions  in  the  form 


atrz  i26' 

Oi  =  1 1 (r j  -  t2(r)cos-j—  . 

The  direction  of  the  bifurcation  is  fixed  now  by  Kj.  Assuming  Kj  *  0,  A(e)  -  ufll  = 

■»  dA  d“A 

Ki and  the  bifurcation  there  is  a  pitchfork  ( - =0and— ■? — ?|)  ate  =  0).  The 

de  d2e 

value  of  K,  is  determined  by  the  soivabilitv  of  (  20c).  and  the  formula  corresponding  to  (23) 

is. 


k. 


4\ 


;2s2<r:uiOo“Oi  +  sKr:|i)On  ) 
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(27) 


and  this  does  not  generally  vanish. 

The  solutions  bifurcating  at  pL01  are  wavetrains  with  wavelength  L  in  an  axially 
infinite  region.  By  developing  the  senes  solution  in  e.  a  finite  amplitude  penodic  wavetrain 
may  be  constructed. 

When  L  — *  <».  fit)l  =  -r  OIL’2),  and  the  solution  of  the  inhomogeneous 
ordinary  differential  equation  for  far;  is  of  OfL2).  As  a  consequence,  the  senes  1 18a) 
becomes  disordered  when  eL~  =  0(1).  The  way  to  deal  with  this  non-uniform  behavior  for 
iong-waves  by  the  metnod  of  multiple  scales  (or  equivalent  methods;  is  weil-known.  In 
the  context  of  the  approach  taken  here,  die  expansion  is  centered  about  the  columnar 
bifurcation  point.  fiu0.  The  procedure,  sketched  in  Leibovich  1 1987),  goes  as  follows. 
Letting  Z  be  the  slow  scale,  with  Z  =  zve.  =  A(Z)<t>o(r)  (at  lowest  order),  then  Opir)  is  the 
principal  eigenfuncnon  corresponding  to  the  eigenvalue  ^K>,  and  A  sausfies  the  equation 


d-A 

JZ- 


«A- 


|iKnA  =  0 


(28a) 


where  *c,  is  defined  by  i  1 8b).  as  before,  but  is  no  longer  restncted  by  the  solvability 
condition  (23).  and 


a  = 


28bi 


^OO')  p  _  '  PlOo") 
,Oo:)  'Of)2} 


The  analysis,  aithougn  bv  a  different  route,  is  essentially  that  of  Beniamin  1 1967  ). 
Equation  (28a)  has  the  solitary  wave  solution 


A  =  a  seen2  [  \'aa/6Zj  =  asech2f  A  aeoc/6  z]  ( 29a) 

provided  aa  >  0  (ana  has  no  such  solution  otherwise).  The  constant  a  in  (29)  is  related  to 
the  parameters  in  (28a):  with  the  level  of  the  extreme  value  of  A  set  to  be  a  and  placed  at  z 
=  0.  as  has  been  done  in  (29).  the  amplitude  a.  or  more  precisely,  ea.  is  linearly  related  to 
the  swirl  rate  by 


la 

-  Max)  -  ~ea- 

a  (3 


(29b) 


In  this  discussion,  it  has  been  assumed  that  £  is  positive,  so  that  waves  of  elevation  or 
depression  depend  on  whether  a  is  positive  or  negative:  since  we  must  have  aa  >  0.  the 
question  devolves  to  the  sign  of  a.  which  is  a  functional  of  the  specifying  flow. 
Furthermore,  since  (3  >  0.  (29b)  shows  that  solitary  waves  may  exist  only  for  A  <  p.,*).  In 
this  parameter  regime,  no  standing  waves  are  possible  (since  they  all  branch  from 
eigenvalues  greater  than  p(K)),  and  by  definition  this  is  a  supercritical  regime.  Thus  (29b) 
makes  clear  that  weakly  nonlinear  solitary  waves  form  only  on  supercritical  flows.  Further 
discussion  of  the  criticality  classification  is  given  in  §5.1. 

The  weakly-noniinear  solitary  wave  solutions  found  by  Leibovich  ( 1970)  from  the 
time-dependent  Korteweg-deVries  equation  are  equivalent  to  those  given  above.  This 
.iltemative  form  effectively  derives  from  the  alternative  parameterization  of  the  time- 
independent  problem  by  wave  speed  (instead of  swiri  level)  mennoned  in  the  Introduction 
and  is 


<t>  =  c a  On(r)  seen2 1 
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w  here  a  is  an  arbitrary  amplitude.  On  is  the  eigenfunction  of  the  linearized  problem 
equivalent  to  (27).  the  c.  are  constants  depending  on  the  base  flow  and  On.  and  ice 
equivalent  to  a  and  (3  in  (28).  For  a  given  base  flow,  the  vaiue  of  ea  determines  the  wave 
amplitude  and  therefore  its  velocity:  or.  alternatively,  the  change  in  the  base  flow  axial 
component  that  wouid  make  the  solution  stationary  .  Changing  /.  (that  is.  the  azimuthal 
component)  instead  will  produce  the  same  resuits  provided  the  swiri  ratio  is  the  same:  so 


'eiecung  an  eigenvalue  and  setting  a.  tor  A)  completely  determines  tms  approximate 
'Olution. 

The  numerical  construction  or  solutions  along  a  branch  can  oe  done  by  a  simple 
conunuation  method,  starting  with  solutions  generated  numerically  by  the  perturbation 
procedure  described  above. 

3.3.  Branch  continuation 

To  "continue  (Kubicek  &  Marek.  1983.  gives  a  good  summary  oi"  continuation  and 
bifurcation  methods)  a  solution  (<DQ,A0)  known  at  a  given  value  of  A  =  A  ,  to  a 

neighboring  value  differing  from  A  incrementally,  we  could  proceed  by  taking  as  an 

initial  guess  for  ct>(A0  +  AA)  in  a  suitable  iterative  procedure,  such  as  Neu'ton  s  method. 


Vl  =  •  L''(®n,A)N«l>„.A)  ,31) 

using  an  appropriate  discretization  for  L  and  N.  TTiis  will  give  a  solution  if  L  is  not 
singular,  that  is.  if  no  bifurcation  or  turning  point  is  encountered,  and  if  AA  is  sufficiently 
small.  To  increase  the  size  of  AA  while  still  providing  a  good  guess  for  the  iteration,  we 
can  proceed  in  the  following  standard  way  by  differentiating  along  solution  arcs  (O(A).A): 


N(<t>(A),A)  =  0. 


dN  GN|  GO  ij  N  | 

— <4>.A)  = — j  - - — !  =0. 

dA  d<l>|0:A  dA  dA  OtA 


Then,  at  any  point  (O.A).  the  "slope  of  the  solution  curve 


•  oO  /GN|  -1  , 

O  - - =  - -  j  PfO)  =  -  L'MO.AiPfO) 

dA  d<D  0;A  j 


(32b) 


assuming  the  linearized  operator  L(d>.A)  is  invertible  with  inverse  L'^O.A).  We  can 
think  of  (31 )  as  a  differential  equanon  for  as  a  function  of  A.  We  use  (31)  and  (32b).  in 
a  numerical  algorithm  described  in  the  appendix,  in  a  predictor-corrector  mode.  We  first 
integrate  (32).  using  a  Runge-Kutta  method  to  arrive  at  an  estimate  for  d>  at  An  -t-  AA 
assuming  a  solution  is  known  at  Ap.  and  then  refine  this  estimate  using  Newton  iteration 
31).  We  begin  each  solution  branch  by  using  a  three-term  perturbation  expansion,  also 
described  in  the  appendix,  which  is  the  discrete  version  of  the  anaivsis  given  in  ij3.2.  This 


procedure  wiil  fail  if  secondary'  Bifurcation  poinis  or  turning  points  are  encountereu  as  a 
given  branch  is  traced,  and  then  a  more  invoived  proceaure.  such  as  tnat  devised  by  Keller 
1 19771  (see  also  Kubicek  &.  Marek.  19831  wiil  be  required.  We  did  not  encounter  such 
complications  in  the  course  or'  our  investigations. 

To  do  detailed  calculations,  we  must  seiect  a  particular  specifying  flow.  We  wiil 
explore  tne  possible  branches  of  solutions  stemming  from  the  following  simple  columnar 
vortex,  which  has  been  previously  treated  by  Leibovicn  ( 1970k 


Wfr)  =  1 


V(r)  = 


exp(-ar^i). 
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This  example  is  known  as  the  Burgers-Rott  vonex.  It  corresponds  to  the  following 
specificauons  of  the  funeuons  ansing  in  $2: 


TTr)  =^r-.  R(y/)  =  2Vy.  ^V)  =  !  -  exp)  -2a\i/) 


->  Ip  1  -  exDi-2aun 

A(y)  =  0;  BCvy.r-)  =  ivy  -  -rr“)  2a  expt  -2au/){ - — - - 1 


PCO)  =  Bf  +  O.r2)  -  B(  ^-".r2)  ;  Q(d>)  =  0. 


3.4  Numerical  Implementation 


<34) 


The  problem  stated  in  (16)  is  discretized  using  central  differences  on  a  rectangular 
mesh  in  the  meridional  plane  tr.z).  Let  <b,  D.  £2  be  the  linite-dimensionai  counterpans  of 
(b.  D-.  Q.  as  derined  in  the  appendix:  1 16)  then  corresponds  to  the  matrix  equauon — 


DO  +  £2(0,  A)  =  0  (35) 

Equation  (35)  is  the  basis  for  the  numencai  treatment.  We  do  not  discretize  steps  of 
the  analytical  procedure  separately:  rather,  we  provide  an  equivalent  analysis  for  the  ap¬ 
proximate  equation  (35).  A  separate  discretization  of  (22).  for  example,  leads  to 
eigenvalues  that  are  slightly  different  than  the  bifurcation  points  of  (35).  and  this  is  enough 
to  prevent  convergence  to  a  solution  branch  in  some  cases.  Even  more  teiling.  if 
eigenvectors  obtained  from  the  algebraic  eigenvalue  problem  are  used  in  conjunction  with  a 
semianalytic  enforcement  of  an  orthogonality  condition  i  for  example,  by  means  of  a 
numencai  quadrature),  then  the  resuit  wiil  not  oe  precisely  onhogonai  in  the  algebraic 
problem,  and  if  the  next  stage  of  the  problem  is  solved  algebraically,  errors  are  introduced. 
We  theretore  re-denve  equauons  ( l~)-(29)  in  the  appendix  for  the  algebraic  system  (35). 
This  ensures  consistency  of  numencai  values  throuenout  the  analysis. 


The  strategy  is  the  same  as  tnat  described  in  §§3.2  ana  3.3.  First,  a  solution  point 
on  a  non-tnviai  brancn  is  sought  using  a  perturoation  expansion  (or.  for  the  soiitary  wave 
branch.  the  weakly-noniinear  solution  may  be  used).  Numerical  integration  or  the  discrete 
anaiog  or  (32b)  continues  the  brancn  away  from  the  bifurcation  point,  and  Newton  s 
iterauons  serve  as  corrector  steps  at  selected  points  along  the  branch. 


4.  Columnar-columnar  bifurcations  and  continuation 

In  this  section,  we  discuss  some  general  questions  about  bifurcations  of  the 
specifying  columnar  flow  to  other  columnar  flows,  and  then  give  numerical  results  for  the 
example  specified  in  (34). 

4.1.  Transfer  of  B-criticality  condition 

Benjamin  <  1962)  has  provided  a  simple  test  to  determine  whether  a  given  columnar 
vortex  is  subcriticai  or  supercritical.  As  Leibovicn  ( 1979)  has  shown,  this  turns  out  to  be 
an  appropriate  test  (for  axisymmetnc  disturbances)  even  though  the  crucial  quantity 
determining  whether  upstream  propagation  of  disturbances  is  possible  is  the  group,  not  the 
phase,  speed.  Subcnncal  flows  can  be  expected  to  be  influenced  by  small  downstream 
disturbances.  This  might  be  true  even  in  flows  that  are  supercritical  according  to  this 
classification  scheme,  since  it  does  not  cover  nonaxiallv  symmetric  perturbations,  but  the 
propagation  characteristics  of  nonaxiallv  symmetric  waves  (see  Leibovicn  et  al..  1986)  is 
more  difficult  to  determine.  Benjamin  s  cridcaiitv  classification  is  important  because  n 
seems  to  be  useful  in  correlating  vortex  breakdown  data  iLeibovich.  1983.  1984),  as 
Squire  ( 1960)  and  Benjamin  i  1962)  had  proposed.  In  particular,  the  evidence  (see 
Leibovich,  1978.  or  L)  indicates  that  flows  upstream  of  vortex  breakdowns  are 
supercritical,  while  the  (mean)  flows  downstream  are  subcriticai. 

To  determine  the  cnncaiitv  condition  of  a  given  columnar  flow,  we  ask  whether  it 
can  sustain  infinitesimal  waves  of  the  form 

o  =x(r)elkz.  (36) 

which  means  that  equation  (16a)  has  a  solution  in  the  form 

<d  +  eo 

tor  infinitesimal  e.  This  leads  to  a  problem  similar  to  that  in  §3.1.  except  that  we  wash  to 
consider  columnar  flows  other  than  the  specifying  flow  ( which  has  =  0),  and.  rather  than 
fixing  the  wavelength  (=  L)  and  searching  for  values  of  A  for  which  the  iineanzed  problem 
is  solvable,  the  question  is  turned  around:  A  is  fixed,  and  we  ask  if  there  is  anv  real  value 
t  k  (=  2~/L)  lor  w  hich  the  linearized  problem  is  solvable.  If  so.  a  standing  wave  w  ith 
wavenumber  k  determined  is  possible  and  the  flow  is  subcnncai.  If  not.  the  flow  is 
supercritical.  Let 


a  .  o~{  )  i  dl  )  ,  aP  , 

LtO.A;  (  i  - - - —  ! A - - d>.r) 
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a<t> 


i  ) 


then  we  may  wnte  the  differential  equation  for  the  small  wavy  pemiroations  as 


d~0 

r)z~ 


k-O  =  LlO.A  )  0. 
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We  now  regard  this  as  an  eigenvalue  problem  tor  k-.  with  A  fixed.  It  is.  due  to 
1 36),  an  ordinary  differential  equation  for  y  in  standard  Sturm-Liouville  form  and  subject  to 
the  boundary'  conditions  y  (a )  =  y{  1 )  ~  0. 

We  know  that  the  specifying  flow  is  supercritical  for  A  <  Mito  (because  there  are  no 
eigenvalues  of  (17)  in  this  range  of  A.  either  corresponding  to  columnar  flows  or  to  wavy 
tlows.  according  to  si  3.1),  ana  subcnttcal  for  A  >  ll,M)  .because,  according  to  the 
observation  of  §3. 1 .  there  is  a  wave  with  some  wavenumber  I  or  each  value  of  A  >  ll,^, 
with  waves  with  indefinitely  long  wavelengths  branching  off  indefinitely  close  to  u<x)-  It 
A  =  pq*),  it  is  clear  that  the  eigenvalue  k-  of  1 38)  vanishes,  and  we  may  summarize  by 
observing  that  the  eigenvalue  k-  of  (  38)  is  negative  if  A  <  zero  for  A  =  p()0,  and 
positive  for  A  >  LnK),  and  that 

dk^ 

— -  >  0  at  A  =  p,K). 
dA 


We  wish  to  characterize  the  criticality  condition  or  the  columnar  vortex  branching 
•  >ff  from  the  specifying  flow  at  A  =  u.K),  on  either  side  or  the  birurcation  point .  I  o 
explore  this,  we  differentiate  (AS)  with  respect  to  A.  to  arrive  at  cm  inhomogeneous 
dy 

equation  tor  — - 
OA 


<ry 

dA 


dA 


dk£  _  dL(O.A) 
dA  ^  dA 
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where. 


dOO.A) 

dA 


o<t> 

oA 
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Since  y  satisfies  i3Si.  (39)  is  solvable  oniy  if 


:  40) 


JL(cD.A) 

/.  — —  X  > 
uA 


We  are  interested  in  (40)  at  A  =  u,*,,  0  =  0.  The  numerator  or  (40)  is 


.X~(  2s.tr.fi/  0  P,(nl) 


where  we  have  used  the  notation  of (21)  and  (32b).  The  only  difference  in  (40).  wtien 
evaluated  at  the  bifurcation  point  for  either  of  the  intersecting  solution  curves  is  the 

•  cii<2 

direction.  <t>.  of  the  particular  solution  curve  aionn  which -  is  evaluated.  On  the 

dA 

specitvinc  tlow.  - —  =  0.  and  on  tne  second  branch  d>  =  y_  .the  eieenfunction.  apan  from 
oA 

an  arbitrary  multiplicative  constant  which  we  take  to  be  unity.  Thus,  along  the  specifying 
tlow.  (40)  reduces  to 


dk^  _  ,PiO')X~)  (42 

dA  "  \X“) 

To  calculate  the  numerator  or  (40)  for  the  nonmvial  bifurcating  solution,  we 
differentiate  (32a)  twice  with  resDect  to  A.  Tlie  result,  when  evaluated  at  d>  =  o.  A  =  u.  , 


LiO.u, „,,<!>  +  2.n, ir.u,,,,/!)*--*-  2p,(riO  = 


and  solvability  tor  O  requires 


<b“{2s  .ir.fi,, ,ii  <b  *  2p.in)>  = 


Commnimi  (43)  with  ( 4 1  i.  we  find  that  on  tne  bifurcatinu  hrancn. 


p ,  ( r  !<t>  “ 


P,'r»X"' 


•which  is  just  the  negative  ot  (42).  Therefore,  since  k-  increases  tnrougn  zero  as  A 
increases  through  on  the  specifying  tlow  brancn.  is  must  decrease  tnrougn  zero  as  A 
increases  through  p<K)  on  the  principal  branching  solution. 

rhus.  we  have  shown  that  the  criticality  conditions  holding  on  tne  specirying  now 
ana  that  holding  on  the  flow  brancning  transcnticailv  from  it  are  transferreu  at  tne 
bilurcanon  point:  this  will  be  illustrated  in  $4.3.  If  we  hold  A  fixed,  our  problem 
conforms  to  that  considered  by  Beniamin  1 1962).  w  ho  labels  each  columnar  rlow  aistinct 
irom  the  specirying  flow'  -  and  there  may  be  more  than  one  depending  on  tne  value  or  A  - 
as  conjugate  '  to  the  specifying  flow.  Benjamin  shows  that,  if  the  primary'  rlow  is 
-upercnncai.  then  all  conjugate  flows  must  be  subcnticai.  This  is  consistent  with  the  local 
resuits  of  this  secnon. 

4.2.  Stability  is  not  transferred 

Stability  or  the  steady  solutions  treated  here  cannot  be  addressed  using  the  BHE. 
-unce  it  contains  no  dynamics.  Instead,  one  must  return  to  the  Euler  euuauons.  As  a  ruie. 
vortex  flow’s  tend  to  less  stable  to  nonaxisvmmemc  perturbations  than  to  axially  symmetric 
ones.  A  bifurcanon  point  is  usually  associated  with  a  transfer  of  stability,  so  that  in  the 
present  case,  one  would  suppose  that  the  specifying  ilow  is  stable  and  the  principal  brancn 
is  unstable  on  one  side  ot  the  bifurcanon  point  A  =  with  the  reverse  being  the  case  on 
the  other  side.  This  is  assured  if  the  eigenvalues  of  the  temporal  linearized  stability 
problem,  deriving  from  the  Euler  equations,  are  simple,  and  provided  they  move  with  a 
nonzero  speed  across  the  imaginary  axis.  In  our  example  (33).  at  least  when  strictly 
restricted  to  axisvmmetric  disturbances,  the  temporal  eigenvalues  are  simple  at  the 
bifurcation  points  ot  the  BHE  equation,  but  bifurcation  need  not  be  associated  with  loss  of 
stability',  which  implies  that  the  temporal  eigenvalues  in  such  cases  are  confined  to  the 
.maginary  axis.  I  his  may  be  seen  from  die  axially-syrnmetnc  Howard-Gupta  equation 
Howard  A  Gupta.  1962).  which  governs  the  temporal  stability  problem.  This  equation  is 
.uenticai  to  1 1 7a i  provided  only  that  tne  axial  veioeitv  of  the  specifying  tlow.  Wi  n.  be 

co 

replaced  by  \V(n  -  c.  where  c  =  —  is  the  (in  genera*  omplex)  phase  speed,  and  <7  =  ;co  is 

the  temporal  eigenvalue.  In  the  example  vortex  (33).  bifurca  :on  does  not  lead  to  ioss  ot 
-lability.  Here  tne  Rayleigh  discriminant  is  positive  and  there  is  no  axial  shear.  Therefore, 
oy  Rayleigh  s  ( 1 SX2)  stability  criterion.  (33)  is  iineariv  stable  to  axisvmmetric  perturbations 
or  ail  A.  On  the  otner  hand,  the  same  stability  condition  must  hold  for  tne  solution 
"ranching  trom  this  point,  at  least  tor  a  iimited  range  of  A.  This  must  be  so  since  the 
mean  zed  stability  characteristics  there  are  determined  bv  the  unperturbed  tlow  at  the 
'(furcation  point,  w  inch  is  the  same  tor  ail  branches  meeuna  there.  Because  tne  velocity 
monies  on  tne  bifurcating  brancn  deform  continuously  with  A.  there  wiii  be  a  finite  A 
ntervai  over  which  the  Richardson  number  criterion  of  Howard  and  Guota.  w  nich 
Generalizes  Rayleigh's  cnterion  to  admit  axial  shear,  continues  to  guarantee  linear  stability 
o  axiaih --ymmetnc  perturbations.  The  eigenvalues,  a.  of  the  temporal  stability  proniem 


.  an  oe  found  from  me  emenvaiues.  u.  of  the  static  bifurcation  nrooiem  in  the  primary  now 
'Oiectea  for  numencai  treatment  (tor  wnicn  Wiri  =  i  i.  since  it  is  easy  to  snow  mat 

A 

a  =  ik  1 1  r  —  ; 

a 

The  eigenvalues  c  are  simple  and  lie  on  the  imaginary  axis  for  ail  A.  and  the  zero 
eigenvalue  is  assumed  when  A  =  u. 

Szen  ( 19881  has  shown  that,  within  the  confines  of  the  Amoid-Casimir  theory, 
these  resuits  may  be  extended  to  weakiy  nonlinear  stability.  (The  interpretation  of  this 
theory  and  its  significance,  as  applied  to  this  proDiem.  is  a  complicated  matter  that  requires 
and  deserves  further  study. ) 

4.3.  Numerical  results 

Figure  1  is  bifurcation  diagram  for  the  columnar  solunons  branching  from  the 
primary  flow  (33.34)  with  a  =  14.  showing  the  principal  and  the  second  bifurcating 
branches.  (The  locations  of  the  points  of  bifurcation  are  given  in  Table  1  in  §5.)  The 
branches  are  described  by  a  measure  of  the  perturbation  axial  velocity.  We  took  this  to  be 
the  extreme  value  (regardless  of  sign;,  and  on  the  principal  bifurcating  branch,  on  which 
the  perturbation  axial  velocity  is  monotonic  in  r.  this  always  occured  at  the  axis  of  rotation. 
On  the  second  branch,  the  perturbation  axial  velocity  if  not  monotonic.  On  this  branch,  a 
discontinuity  appeared  in  the  bifurcauon  diagram  based  upon  this  measure  described:  for  a 
range  of  swiri  parameter,  the  extreme  value  occured  on  the  axis,  but  shifted  to  a  point  off 
:he  axis,  w  here  the  perturbauon  axial  velocity  was  of  opposite  sign.  This  is  illustrated  by 
:wo  sets  of  points  in  the  figure,  one  as  described,  and  the  other  (smooth)  set  arrived  at  by 
plotting  oniv  the  perturbation  axiai  velocity  on  tne  axis. 

(FIGURE  1  ABOUT  HERE) 

The  principal  columnar  branch  plays  an  important  role,  as  will  be  seen  7 
contrast,  the  physical  significance  of  the  second  and  higher  columnar  branches  u,  ungear  - 
for  the  larger  amplitude  perturbations  on  tnese  branches,  axial  flow  reversals  are 
necessarily  accompanied  by  internal  zeros  of  the  swirl,  and  hence  instability  according  to 
the  Rayleigh-Synge  (  1933)  criterion  (abbreviated  subsequently  by  R-S).  This  may  be  seen 
m  figure  2.  which  shows  profiles  from  three  points  on  the  second  branch  of  figure  i.  Here 
the  bifurcation  point  occurs  at  A  =  quo.  and  the  profiles  are  drawn  for  increasing  vaiues  of 
V  Ltii).  One  can  aiso  see  from  tne  protiles  for  axial  velocity  how  the  discontinuity  in  the 
bifurcation  diagram  discussed  in  the  previous  paragraph  arises.  We  wail  not  further  discuss 
columnar  branches  other  than  the  principal  branch. 


-:o- 


(FIGURE  2  ABOUT  MERE) 


For  A  <  the  pnncipai  branch  shows  a  developing  wake-like  axiai  velocity 
profile  as  the  swiri  parameter  A  is  decreased  from  the  brancn  point  u,v>  as  may  be  seen  in 
figure  3.  The  swiri  velocity  is  distorted  as  weil.  with  the  peak  swirl  moving  outwards 
relative  to  that  in  the  primary  vortex.  Both  of  these  characteristics  are  qualitatively  like  the 
lime-averaged  profiles  measured  by  Garg  and  Leibovich  ( 1979)  (further  analysis  of  this 
data  is  given  by  L  and  by  Leibovich.  1978)  downstream  of  vortex  breakdowns,  which  are. 
like  the  solutions  here,  wake-like  and  subcritical.  If  the  A/pgo  is  decreased  below  a  value 
of  about  0.5,  the  swirl  velocity  develops  an  internal  zero,  and  the  branch  will  become 
unstable  according  to  the  R-S  criterion.  We  have  nevertheless  continued  to 

(FIGURE  3  ABOUT  HERE) 

follow  the  branch,  being  curious  to  know  if  it  couid  be  continued  to  zero  swiri.  or  whether 
it  would  turn  around.  We  found  that  it  can  be  conunued  to  zero  swirl,  and  the  curvature  of 
the  magnified  pan  of  the  bifurcauon  diagram  of  figure  1  indicates  that  the  axiai  speeds  get 
large  as  A/p^g  —*  0.  The  axial  profiles  shown  for  A/p^g  =  0.001  suggest  that  the  limit  flow 
becomes  disconunuous.  with  a  vortex  sheet  forming  in  the  intenor.  The  resoluuon  of  the 
singular  limit  behavior  as  A/p^g  — ■>  0  requires  special  treatment,  and  §4.4  is  devoted  to  this 
question. 

On  the  A/q+x)  >  1  side  of  the  principal  branch,  the  axial  velocity  profiles  are  jet-like, 
and  the  peak  swiri  moves  towards  the  axis  relative  to  that  in  the  primary  vortex.  Examples 
are  shown  in  figure  4.  No  tendency  towards  R-S  instability  occurs  on  this  side  of  the 
bifurcation  point.  These 


(FIGURE  4  ABOUT  HERE) 

velocity  profiles  not  only  resemble  the  profiles  measured  by  Faler  and  Leibovich  (1977.8) 
and  in  the  references  cited  in  the  previous  paragraph  for  flows  well  upstream  of  vortex 
breakdown,  they  aiso  can  be  accurately  fitted,  as  can  the  wake-like  soiuuons  previously 
discussed,  by  the  same  exponential  functions  used  in  those  references.  We  note  further 
that  the  experimental  data  shows  the  upstream  flow  to  be  not  only  jet-like,  but  supercritical. 
Thus,  the  primary  vortex  with  uniform  axiai  velocity  generates,  through  its  pnncipai 
branch,  vortices  of  the  same  character  as  those  found  on  both  upstream  and  downstream 
sides  of  experimentally  observed  vortex  breakdowns,  so  far  as  the  shapes  of  the  profiles 
and  their  criticality  conditions  are  concerned. 

For  A  <  P(M),  the  Dnmarv  vortex  is  supercritical  and  the  pnncipai  branch  is 
subcritical.  and  according  to  the  general  theory  of  §4. 1.  these  charactenstics  should  be 
exchanged  w  hen  A  >  u,H).  We  have  tested  this  by  computing  the  generalized  Froude 
number.  N,  proposed  by  Beniamin  ( 1962).  This  is  defined  to  be 


N 


'.vhere  c*  ana  c.  are.  respectively,  the  maximum  and  minimum  phase  speeds  of 
infinitesimal  axisymmetnc  waves  of  extreme  iength  propagating  on  the  vortex.  If  N  >  1 ,  a 
vortex  is  supercritical,  and  if  N  <  1  it  is  subcritical.  We  have  computed  N  on  both  the 
pnmarv  vonex  and  the  pnnctpal  branch,  and  the  resuits,  given  in  figure  5.  confirm  the 
general  theory  on  exchange  of  criticality  at  the  bifurcation  point. 

(FIGURE  5  ABOUT  HERE) 


4.4.  Singular  limit  of  zero  swirl 

For  the  specitving  rlow  (34)  that  we  have  been  using  as  example,  the  zero  swiri 
limit  A=0  is  irrotational.  and  the  constant  speed  flow  u  =  r  r;  is  unique.  Therefore  no 

differentiable  soiuuons  exist  except  for  the  specifying  flow.  On  the  other  hand,  our 
numerical  results  suggest  that  the  limit  A— » 0+  along  the  principal  columnar  branch 
develops  a  strong  shear  layer  tending  in  the  limit  to  a  vonex  sheet  separating  ru4&f«rewise 
irrotanonal  limiting  flows,  in  each  section  of  which  the  flow  is  oppositely  directef^iil  of 
uniform  infinite  speed!  We  explore  this  bizarre  situadon  by  an  asymptotic  deve’ojjaiilit  that 
seems  to  fit  with  the  numerical  findings,  and  which  gives  the  asymptotic  dependcdfceof  the 
velocity  levels,  and  of  the  shear  layer  location  and  thickness  on  A. 

It  is  convenient  here  to  work  with  the  total  streamfunction  vir.  governed  by  <11 ). 
rather  than  the  perturbation  from  the  specifying  flow.  For  the  specifying  flow  1 34),  the 
.unit  or  equauon  (11)  when  A=0  is  D-ur  =  0.  with  a  solution  that  is  iinear  in  r~.  We  assume 
the  existence  oi  a  discontinuity  in  this  iimit.  located  at  r=r~.  at  which  the  stream! unction 

reaches  its  minimum  value,  wr*)  =  -  ':/*(  A) .  The  analysis  of  this  section  in  simplified  bv 
the  change  of  variable 


n  =  T  r*-. 


Then  we  may  write  (using  the  condition  that  ti/  =  0  at  rj  =  0  and  \j  =  ~  at  p  =  and 
defining  ty*  =  !n-i 


'./  =  ('  p.  tor  q  -  q«  <  i) 

';/  =  D  ( q  -  -b  -  - .  for  q  -  rj «  >  i ). 


1 43) 


For  A  very  small  but  not  zero,  we  assume  the  existence  of  a  single  internal  layer, 
centered  on  p*  and  with  a  thickness  6(A)  that  tends  to  zero  as  A— >  0+.  that  toms  these  two 
constant  speed  solutions.  Furthermore,  the  numerical  resuits  suggest  that  xj*( A)  — *  «=  ;is 
A— >  0.  and  we  also  assume  this. 

The  solution  to  the  outer  problem  nas  been  descrmed  above,  and  now  we  seek  a 
structure  to  the  internal  boundary  layer  separating  the  two  irrotauona]  regions,  it  is 
important  (and  easy  to  show)  that  D2g/  =  0  to  ail  algebraic  orders  in  the  small  parameter  in 
the  outer  regions,  so  that  the  full  outer  expansions  retain  the  form  (45),  with  the 
coefficients  of  r)  being  functions  of  A.The  point  p*.  and  the  constants  C  and  D  occunng  in 

the  outer  problem  are  not  yet  known  and  must  be  determined  by  matching  with  this  internal 
boundary  layer.  The  argument  is  reminiscent  of  activation  energy  asymptotics  tsee 
Buckmaster  and  Ludford.  1982).  Stretch  the  radial  scale  near  p  =  p*  .  by  taking 

p  =  p*  -r  OX  =  p*(  1  +  uX  )  .  (I  =  — -  ,4f>) 

T”|  * 


where  6  is  the  length  scale  appropriate  in  the  layer  and  it  is  assumed  that  6(A)  *s 
A— >0.  Near  p*  the  streamfunction  is  continuous  and  the  appropriate  scale  for  it  is- jL(A), 
<o  we  wnte 

\|/(p*  +  6X)  =  -  vg*(A)[l  +  e(A)  y(X;A)]  Jr  <47-) 

in  the  layer,  where  the  asymptotically  small  parameter  £(  A),  like  the  parameters  6(A).  and 
>4*(  A),  remain  to  be  identified.  From  the  definition  of  -  './*(  A)  as  the  extreme  vaiue  of  vi/. 
we  must  have 


y(0;A )  =  y  (0:A)  =  0  (48) 

dv 

where  (  f  =  — .  and  we  also  have  y  <  0  for  ail  X  .  We  substitute  the  ansatz  (45-47)  into 
(11)  and  invoke  ( 34).  This  yields 

2aA  6“ 

V  +  -  -  <r{\)  =  0 

t'U/* 


where 


O'fv)  =  (— - 


-«W,U  +fVI,  .  --UlK.ll  rtVI 

e  (  1  -  >: 


_n.  _u*t  i  -evi 

The  distinguished  limit  arises  tor  2a\i/*  £  ana  2a A6~e"+au/*  '2r\-*  doth  0(  1 ).  ana  we 
therefore  set 


2a\i/*  £ 


l  7  A  51  ^  4(11(1/*  i 

=  da-Ao^-e  /2q*=l. 


which  leads  to  the  equation 


(49) 


2r|*ot£ 

4- - 

1+  £v 


(50) 


with  exponentially  small  error.  We  may  now  iook  tor  a  solution  to  this  inner  Droblem  in 
the  form  of  power  senes  in  the  small  parameter.  It  is  more  convenient  in  the  analysis  to 
regard  £  as  the  controlling  smail  parameter  rather  than  A.  and  so  we  take 


V  =  yo+  evi  +  ... 

(51) 

U  =  Ui£  4-  M.2£^  4-  ... 

(52) 

~~  =  y0  +Yi£  +  ••• 

(53) 

and  the  slightly  unconvennonal  form  of  the  last  expansion  makes  the  matching  with  the 
■  niter  solution  simpler. 

The  equations  for  first  two  coefficients,  yo  and  vj.  of  (5 1 )  are 


y  o"  + 


(54) 


yf  + 


1 55) 


The  solution  to  (54)  satisfying  t4N)  is 


yt)(X)  =  in  scch(X). 

and  the  solution  to  (55)  is 


56) 


. '  i- 


yi<X)  =t[(J-i(X  -  tann  X)  -  —  Xtann  xj. 

—  '  / 

•  n 

We  now  match  the  streamfunction  in  the  inner  region  with  the  two  outer  regions. 
Matching  to  second  order  in  e  yteids 


— 1 —  =  i2i-E(l+ai)  +  v,t;X  1 58a  i 

2p* 

C  =  - — - — H+eln2--U:>  <  58b) 

2aep* 

D  = - - — <-l  +  e  { 1- in  2)+  e^b->)  (58ci 

2a£T|* 

u  =  e  +  e^  (  7-f  tcx  -  in  2)  <58d) 

_  -+■ 

5  =  ri*(a,  ( 58e ) 

and  from  these  we  can  find 

A  =  2p*(2ot8  expt2ode))’“.  l58f) 


The  numbers  y->  in  (58a)  and  b?  in  (58c)  are  undetermined  at  this  order,  and  therefore 
our  solution  is  completely  determined  only  to  one  order  less  than  we  have  shown  -  on  the 
other  hand,  it  is  necessary  to  match  at  the  level  shown  to  accurately  determine  the  solution 
:o  that  order  ti.e..  to  within  an  error  of  Ole)). 

The  composite  expansion  for  \j  constructed  from  (45-47.56  58>  is 


1  r  p 

vj  = - \  (i  +  t  in  2)  -4-  Hip*  -  n) 

2ae  p* 


(59) 


n  ‘  T 


(-1  +  e  ( 1-  In  2-r- 


p* 


■  flip  -  p*l 


,  i  .  n  -  n  * 

e  lnrsecn  ( - )  -r 


p  -  n* 


! 


t  Hip  -  p  < );  v/(r*i  = - 

2ae 


a  here  Hit)  is  the  unit  function.  Hit)  =1  for  i  >  1 1.  1 1(0)  =  -.  ana  Hit)  =  o  for  t  •  t>. 


Equations  (58 1  give  6.  r*,  and  A  as  a  funcnon  or  the  parameter  e.  and  from  tnem 
the  dependence  or  6  ana  r*  on  A  can  be  evaluated  to  vieid  the  desired  asymptotic 
relationship  between  these  parameters  as  A— *0.  As  seen  in  figure  o.  the  differences 
between  the  computed  and  the  asymptotic  values  are  only  a  few  percent  for  A  =0.01  The 
asymptotic  solution  (59)  is  compared  in  figure  7  to  the  computed  streamiunction  and  axial 
velocity  for  A=  1.769- 1  ()"■*.  The  agreement  is  quite  good,  considering  tne  fact  that  vi/*= 
0.13.  not  yet  large  as  the  analysis  assumes,  meaning  that  the  asymptouc  relation  is  good  far 
beyond  its  expected  region  of  iv|/*l»  i . 

(FIGURES  6  &  7  ABOUT  HERE) 

The  azimuthal  velocity  component,  v,  is  exponentially  small  except  in  the  shear 
layer.  The  shear  layer  is  a  concentrated  region  of  voracity,  with  both  axial  voracity  arising 
from  the  swiri  as  weil  as  azimuthal  voracity  arising  from  the  variation  in  axial  velocity 

across  the  layer  present.  As  A-^().  the  swirl  component  v  —  0  in  the  layer  as  \  A  /lln  A|. 
which  is  very  much  slower  than  A. 

5.  Periodic  wavetrains 

The  specifying  flow  is  subcntical  for  A  >  for  any  A  in  this  range,  an 
infinitesimal  standing  wave  is  possible,  with  wavelength  depending  on  A.  Fixing  the 
wavelength  at  L.  we  continue  infinitesimal  standing  waves  to  finite  amplitude  for  the 
specifying  flow  ( 34)  with  a  =  2  and  a.  =  14  using  the  methods  previously  described.  The 
larger  value  of  a  was  used  in  an  earlier  study  bv  Leibovich  (1970),  because  it  provided  a 
good  fit  to  Harvey's  ( 1962)  experimental  data.  The  eigenvalues  ppm,  determined 
numerically  as  described  in  the  Appendix,  are  given  in  Table  1  for  a  =  2  for  L  =  o  and  20. 
and  for  a  =  14  for  L  =  a.20.100. 


« 

L 

Poo 

dot 

Pi  0 

- 

6 

1.8433 

1.9693 

0.2743 

20 

1.8433 

1.8547 

0.2743 

14 

0 

.17390 

.17986 

.01147 

20 

.17390 

.17445 

.01147 

100 

.  17390 

.17392 

.01147 

Table  i .  Eigenvalues  ot  the  linearized  problem  for  selected  values  of  L  anu  a. 
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Figure  »  is  the  bifurcation  diagram  for  the  case  oc  =  2  ana  L  =  iO.  The  diagram  for 
a  =  14  is  qualitatively  similar,  but  the  separanon  between  the  curves  for  the  columnar 
branch  and  the  first  wavy  branch  cioses  too  rapidly  to  oe  conveniently  illustrated  in  a 
drawing.  The  soiid  line  gives  the  coiumnar  (00)  brancn.  open  circles  the  fundamental 
wavy  branch  (01).  and  the  close  circles  the  first  harmonic  wavy  branch  (02). 

(FIGURE  8  ABOUT  HERE) 

We  next  explore  the  waveform  for  three  flows  on  the  fundamental  wavy  branch  for 
the  case  a  =  14.  L=10.  The  streamfunction  at  a  fixed  value  of  r  =0.25  is  plotted  in  figure  9 
as  a  function  of  z  over  one  wavelength  .  for  three  values  of  A/Uoo-  A/u<x)  increases,  the 
wave  trough  becomes  increasingly  sharp  and  concentrated,  and  the  wave  crests 
increasingly  broad  and  flat.  The  same  trend  is  seen  if  the  fundamental  wavy  branches  for  a 
sequence  of  flows  of  increasing  wavelength  L  are  sampled  at  fixed  A/ mo- 

( FIGURE  9  ABOUT  HERE) 

Figure  10  shows  streamlines  projected  onto  a  meridian  plane  for  increasing  values 
of  A/|ioo  for  the  case  L  =  b  and  a  =  1 4.  The  deceleration  of  the  upstream  flow  caused  by 
the  wave  is  apparent  at  the  smallest  value  of  A/|loo  shown.  The  other  two  streamline  tields 
reveal  a  region  of  closed  streamlines,  with  the  size  of  the  recirculation  region  growing  with 
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(FIGURE  10  ABOUT  HERE) 


The  changes  in  waveforms  as  either  (A  -  p,:)(pp<io  or  L  increases  are  illustrated  in 
figure  9.  Increases  in  either  of  these  parameters  appear  to  produce  a  wave  shape  with  a 
very  sharp  trough,  in  w'hich  there  are  strong  axial  accelerations,  rapidly  tending  to  a  broad 
flat  crest.  Over  most  of  its  extent,  this  broad  crest  is  an  essentially  coiumnar  flow,  but 
distorted  considerably  from  the  primary  columnar  flow.  These  features  are  characteristic  of 
a  solitary'  wave,  with  L  =  =».  on  a  columnar  flow  different  from  the  specifying  flow. 

Velocity  profiles  at  the  wave  trough  are  given  in  figure  1  la.b  for  a  two  swiri  levels 
large  enough  to  cause  a  region  of  reversed  axial  flow  to  appear.  Figure  1 1c  shows  the 
difference  between  the  axial  velocity  at  the  trough  and  the  axial  velocity  at  the  crest,  where 
the  flow  is  neariv  columnar  and.  as  wail  be  shown,  virtually  indistinguishable  from  the 
principal  conjugate  flow.  This  is  a  perturbation  caused  to  the  principle  coniuuate  flow  by 
the  disturbance  concentrated  near  tne  wave  trough,  rather  than  a  perturbation  to  the 
<pecitying  flow.  The  profiles  in  i2c  are  similar  to  the  perturbation  axiai  velocity  in  solitary 
waves  on  the  snecifvinsi  flow. 


(FIGURE  !  i  ABOUT  HERE) 


The  observations  we  have  made  about  the  apparent  appearance  ot  solitary  wave 
behavior  will  now  be  put  to  quantitative  tests. 

The  characteristic  ienath  ot  a  solitary  wave  can  be  measured  by,  say,  its  half-height 
length.  According  to  the  weakly  noniinear  soiiton  solution  (29),  the  half-height  length  (or 
anv  other  measure  of  the  soiitarv  w-ave  length)  scaies  with  the  inverse  square  root  ot  the 

w'ave  amplitude,  l/y/ea.  or  alternatively,  with  1/V I  A— |dooh  at  least  tor  IA  -  Moo* 
sufficiently  small.  The  axial  velocity  perturbation  at  the  origin  is  proportional  to  the  wave 
amplitude,  £a.  In  figure  12  we  have  plotted  the  half-height  calculated  from  our  numerically 
determined  solutions  on  the  first  wavy  branch  for  several  wavelengths  L  =  (6,20,100), 
against  the  axial  velocitv  perturbation  at  the  on  gin  to  test  to  see  whether  the  amplitude- 
length  scaling  appropriate  to  the  weakly  noniinear  solitary  wave  is  approached  by  a  wave  ot 

(he  computed  penodic  wavetrain.  Waves  exhibiting  that  scaling  will  have  a  slope  of  -  -e  on 
this  log-log  plot,  and  the  soiid  line  drawn  has  tnat  slope.  The  wavetrain  with  penod  L  =  ti. 
given  by  the  open  circles,  deviates  substantially  from  the  -  —  slope  both  tor  small  and  large 

wave  amplitude.  The  longer  waves,  however,  accurately  display  the  solitary  wave  scaling 
for  amplitudes  ranging  from  fairly  small  values  to  quite  substannal  ones.  Marked 
deviations  from  the  weakly  nonlinear  solitary  wave  scaling  occur  for  tor  very  small  axial 
velocity  perturbations  and  for  axial  velocity  perturbations  ot  0(1)  and  higher.  As  will  be 
seen  in  the  next  section,  strong  nonlinear  effects  distort  the  solitary  wave  on  the  primary 
tlow  in  the  same  way,  and  begin  to  substantially  modify  the  weakly  nonlinear  scaling  at 
about  the  same  level  of  axial  velocity  perturbation.  These  deviations  from  the  weakly 
nonlinear  scaling  do  not  signal  departures  from  solitary'  wave  behavior,  but  rather 
transitions  to  a  strongly  nonlinear  solitary  wave  regime.  We  note  that  lor  large 
perturbations,  the  L  =  n  case  falls  on  the  same  curve  as  the  longer  waves,  implying  that  the 
penod  of  this  wave  is  not  long  enougn  to  exhibit  solitary'  wave  behavior  at  smail  ampiitude. 
but  that  it  does  develop  strongly  nonlinear  solitary  wave  behavior  at  large  amplitude. 


(FIGURE  12  ABOUT  HERE) 


We  conclude  from  these  considerations  that  one  wavelength  ot  a  penodic  wavetrain 
rapidly  approaches  a  soiitarv  wave  as  the  wave  ampiitude  increases  above  a  modest  level. 
The  resulting  motion  may  be  charactenzed  either  as  weakly  or  as  strongly  nonlinear  solitary 
waves,  depending  upon  amDiitude.  Funhermore.  the  columnar  tlow  to  which  these 
solitary  waves  tend,  at  distances  from  me  station  of  maximum  amplitude  large  compared  to 
(he  half-length  Lip,  is  not  the  specifying  tlow-.  but  the  principal  conjugate  tlow.  This  point 
is  illustrated  in  figure  l.\  Figure  ibu  shows  tne  difference  between  the  streamrunction  ot 
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principal  conjugate  flow  and  that  ot  the  first  wavy  nrancn  as  a  function  ot  r  tor  waveiengtn 
L  =  6  for  three  vaiues  t().()5.  0.25.  ana  1.0)  of  (A  -  corresponding  tne  waves  of 

increasing  amplitude.  At  the  largest  t  A  -  u<K)y  p.**,  thence  the  largest  amplitude i.  die 
difference  is  bareiv  detectable.  Figure  13b  shows  the  same  tendency  as  L  increases  with 
the  wave  amplitude  fixed. 


(FIGURE  13  ABOUT  HERE) 


6.  Solitary  waves  on  the  specifying  flow 

When  the  specifying  tlow  is  slightly  supercritical,  a  weakly  nonlinear  solitary  wave 
<29)  is  possible.  A  diagram  summarizing  the  numerical  continuation  of  this  solitary  wave 
solution  brancn  to  more  strongiv  supercritical  conditions  (A  decreasing  from  u,K);  note  that 
the  scaled  distance  from  the  branch  point.  IA  -  u^l/poo,  ranges  from  zero  to  a  maximum  of 

unity;  is  given  in  figure  14.  The  diagram  superposes  two  measures  of  the  wave  disturbance 
of  the  specifying  columnar  tlow.  the  maximum  perturbation  axial  velocity  (wmaxj  at  the 

plane  of  symmetry  z  =  0.  and  the  perturbation  axial  velocity  on  the  axis  at  this  plane 
(w'(0,0)).  The  two  measures  agree  for  values  of  IA  -  pool/p^o  as  large  as  0.8.  For  larger 
values  of  this  parameter,  the  point  at  which  the  perturbation  axial  velocity  is  a  maximum 
lifts  off  of  the  rotation  axis. 


(FIGURE  14  ABOUT  HERE) 

The  shift  of  the  point  of  maximum  axial  velocity  disturbance  may  also  be  seen  in  the 
axial  velocity  profiles  at  the  symmetry'  plane.  Profiles  of  the  complete  axial  velocity 
comDonent  are  drawn  in  figure  15  for  four  vaiues  of  A.  Three  ot  the  profiles  include 
negative  vaiues  of  w.  which  implies  the  existence  of  a  region  of  closed  streamlines 
containing  reversed  axial  velocities.  When  the  maximum  perturbation  lifts  off  the  axis,  a 
high-speed  upstream-directed  iet  forms  in  the  intenor  of  the  recircuiauon  region,  and  the 
dividing  streamline  develops  a  dimple  at  the  axis  and  is  no  longer  convex.  We  are  unawaire 
of  observations  of  such  a  phenomenon,  and  believe  it  to  be  physically  unrealizable. 

(FIGURE  15  ABOUT  HERE) 

Projections  of  the  streamlines  onto  the  mendian  plane  are  shown  in  figure  16. 

These  plots  show  the  emergence  of  the  recirculation  region.  We  have  found  that  the  tlow 
field  is  represented  with  reasonable  accuracy  by  the  weakly  nonlinear  solution  (29)  for 
waves  leading  to  axial  velocity  perturbations  strong  enough  to  cause  tlow  reversal.  This  is 
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a  significant  rinding,  since  we  may  men  capture  the  essenuais  or' our  numerical 
computations  or  a  strongiv  penurbed  flow  with  the  simpie 

i  FIGURE  16  ABOUT  HERE) 


r’ormuia  (29).  To  snow  the  levei  of  agreement,  we  have  compared  (for  a  =  i4)  w'(O.O) 
from  our  numencai  computations  with  the  approximation  (29).  This  is  found  in  figure  17. 
together  with  a  comparison  of  the  dependence  of  wave  half-length  with  wave  amplitude 

(measured  by  w  (0.0)  with  the  -  ^ power  law  dependence  obeyed  by  the  weakly  nonlinear 

solitary  wave.  The  weakly  nonlinear  solution  (29)  overpredicts  the  wave  amplitude  and 
length,  but  the  differences  are  less  than  10%  for  wave  amplitudes  large  enough  to  cause 
stagnation  and  reversed  axial  flow. 


(FIGURE  17  ABOUT  HERE) 

Contours  of  the  perturbation  streamfunctions.  as  predicted  by  numencai  computation 
and  by  the  weakly  nonlinear  approximation  (29),  are  snown  in  figure  i8.  We  judge  the 
agreement  to  be  qualitatively  good  for  ail  three  values  of  1A  -  u^'/flno  shown, 
quantitatively  good  for  I A  -  (%)!/(%)  =  0.1.  and  acceptable  for  some  purposes  for  the 
higher  values  of  IA  -  (J-oo^M-oo-  It  is  worth  noting  that  a  stagnation  point  first  appears  in  the 
flow  for  IA  -  Poo'/hoo  =  0.25,  so  the  three  cases  shown  in  figure  17  range  from  moderately 
to  strongly  nonlinear.  &3K 


7.  Discussion  and  Conclusions 


We  have  shown  here  the  connections  between  fully  nonlinear  standing  periodic 
wavetrains  and  solitary  waves  and  the  underlying  columnar  flows.  From  a  given  primary, 
or  specifying  '.  columnar  flow,  other  columnar  flows,  solitary  waves,  and  periodic 
wavetrains  may  be  constructed.  The  solitary  waves  exist  only  when  the  primary  flow  is 
supercritical,  a  condition  that  arises  when  the  swiri  rate  or  axial  voracity  in  the  primary 
flow  is  less  than  an  easily  determined  critical  value.  Periodic  wavetrains  exist  oniy  when 
the  primary  flow'  is  subcruical.  which  arises  when  the  swirl  rate  exceeds  the  critical  value. 
On  the  other  hana.  the  periodic  wavetrains  rapidly  attain  the  characteristics  of  solitary 
waves  as  the  swiri.  and  with  it  the  wave  amplitude,  increase  for  fixed  wavelength:  or  as  the 
wavelength  increases  at  fixed,  but  finite,  amplitude.  These  solitary  waves  do  not  propagate 
on  the  primary  subcntical  columnar  flow.  Instead,  the  flow  far  from  the  wave  center  is  the 
principal  conjugate  columnar  flow.  The  latter  flow,  in  turn,  connects  to  the  primary  flow  at 
the  critical  swiri  level,  and  as  we  show’  is  supercritical  w  hen  the  primary  flow  is  subcriticai. 
Thus  the  requirement  that  flow  upstream  of  solitary  waves  must  be  supercritical  is 
maintained. 

The  simple,  partly  analytical,  formula  for  weakly  nonlinear  solitary  waves  is  shown 
to  fit  the  numencai  data  for  fuilv  nonlinear  solitary  waves  verv  well  for  a  substantial  ranee 


of  amplitudes.  The  errors  assoeiateu  with  this  lit  are  reiativeiv  small  even  ior  waves  w  un 
amplitudes  large  enough  to  cause  staenation  points  and  reversed  axial  How  to  occur. 

When  stagnation  points  form,  ana  with  them  recirculation  regions  of  closed 
Ntreamiines.  we  must  lace  the  special  quesnons  concerning  the  interpretation  of  the  resuits. 
This  is  due  to  the  well-known  nonuniqueness  or'  steady,  axially  symmetric,  inviscid  flows 
with  closed  streamlines.  When  closed  streamsurfaces  exist,  the  specification  or  the 
voracity  distnbuuon  by  functional  forms  for  H'fy)  and  Frig)  determined  by  the  upstream 
flow  need  not  be  continued  into  the  region  of  closed  streamlines.  In  fact,  if  the  flow  is  to 
be  steady  and  to  be  the  limit  of  a  viscous  flow  as  the  viscosity  vanishes,  the  voracity  in  the 
recirculation  region  must  satisfy  the  constraint  found  by  Prandtl  f  1 904')  and  Batchelor 
( 1956),  W'hich  requires  F(v|/)  =  0  and  H'(\g)  =  constant.  One  might  think  that  a  solution 
with  closed  streamlines,  ignoring  the  Prandtl-Batchelor  tor  PB)  condition,  can  be  made 
consistent  with  it  by  recomputing  the  flow  inside  the  dividing  streamsurface.  using  the  PB 
criterion  to  fix  the  interior  voracity  distribution,  but  maintaining  the  extenor  flow  and  the 
shape  of  the  dividing  streamsurface.  This  generally  cannot  be  done,  however,  while 
balancing  the  intenor  and  extenor  pressures.  Thus,  if  one  insists  on  satisfying  the  PB 
condition,  flows  computed  with  an  arbitrary  specifying  flow  upstream  must  be  discarded  if 
recirculation  appears  unless  a  free  streamline  type  of  problem  is  solved  in  which  the 
boundary  shape  between  the  external  flow  and  the  internal  (PB  flow)  is  pan  of  the  problem 
(see  Leibovich,  1968  for  an  example  of  such  a  construction). 

On  the  other  hand,  the  PB  enterion  fails  if  the  flow  is  not  truly  steady.  ( I*  may 
fail  in  other  circumstances,  as  described  by  Batchelor.  1956.)  It  is  eac"  to  imagine  dearly 
inviscid  flow  developing  due  to  extemai  forcing  of  vanous  kinds,  and  then  settling  1g1>  a 
phase  of  very  slow  change.  If  closed  streamlines  are  present,  then  viscous  effect::  xiws  act 
and  cause  the  flow  to  vary  with  time.  But  this  development  is  very  slow,  and  so  i..#ie  is 
interested  in  time  intervals  shon  compared  to  the  viscous  nme  (of  order  Rr/v.  where  R  is  a 
length  scale  characterizing  the  recirculation  region),  then  the  PB  condition  does  not 
constrain  the  voracity  distnbuuon.  The  PB  criterion  also  fails  if  the  axial  symmetry  is  lost 
'even  if  the  nonaxiailv  symmetnc  component  of  the  flow  is  infinitesimally  weak)  since 
there  is  then  fluid  exchange  across  the  nominal  closed  streamsurface.  Our  interest  in 
vortex  breakdown  leads  us  to  contemplate  flows  in  which  001)1  of  these  conditions  (external 
forcing  dnving  flow  axially-svmmetnc  development  of  d  sed  streamline  regions  on  inertial 
time  scales,  followed  by,  or  coincident  with,  symmetry  breaking  instability)  are  active.  If 
an  axially-svmmetnc  recirculation  zone  existed  and  then  was  broken,  the  fluid  exchange 
across  the  nominal  boundary  would,  in  our  view,  create  an  intenor  voracity  distnbuuon 
that  is  not  inconsistent  with  that  in  the  extemai  flow  upstream.  Thus  the  extemai  form  for 
H'ri|/)  would  be  reasonable  in  the  intenor  (but  the  enure  flow  would  of  course  be  perturbed 
by  the  asymmetric  motions),  and  a  reasonable  form  for  the  circulation  in  the  intenor  would 
be  not  Frig)  but  -  Frig).  This  alteration  is  dynamically  compatible  with  the  flows  computed 
here  tas  pointed  out  in  L)  and  produces  a  flow  with  an  intenor  swirl  sense  in  agreement 
with  that  in  the  extemai  flow,  which  is  certainly  required  if  there  is  exchange  of  fluid. 

Ill  us.  it  is  our  view  that  the  flows  calculated  here  having  recirculation  regions  are  sensible. 
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although  possibly  too  simplistic,  models  of  real  vortex  tlows  with  stagnation  points  ana  a 
>emblance  ot' a  recirculation  region  (albeit  a  broken  one).  We  are  in  the  process  of 
exploring  the  breaking  of  the  symmetry  of  these  flows  produced  in  this  paper,  and  intend  to 
report  on  that  investigation  in  the  future. 
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Appendix:  Numerical  Implementation 


A.l  Discretization 

A  finite-difference  discrenzanon  is  done  on  equauon  1 16).  written  in  the  form: 

,  .  rl  /I  rjO  y  rj~<  t>  I  A  1  | 

D-<3>  +  i2fcD.r:A)  =  r— •'  - -  -  L2(&s:A)  -  t)  A 

orv  r  or  ! 

Die  finite-difference  approximation  corresponding  to  this  form  of  the  radial  derivatives  has 
the  property  that  the  contribution  of  each  cell  boundary  to  the  circulation  around  the  ceil  is 
the  same  for  the  two  cells  adjacent  at  that  boundary:  the  result  is  that  any  group  of  cells 
satisfy  the  Stokes  theorem  when  individual  ceils  do.  similar  to  the  flux  consistency  in 
■'conservative''  discrenzation  of  the  Navier-Stokes  or  the  energy'  equauon. 

We  expect  solutions  that  have  sharp  axial  gradients,  and  non-uniform  grid  spacing 
may  be  necessary'  in  the  z-direction.  This  is  done  by  defining  a  computational  coordinate  C. 
reiated  to  the  physical  coordinate  z  by  z=f(C)  (ffOi^O);  the  z-derivative  in  i  A 1 )  becomes: 

_ i  d  (  i  odd 

o2 Z  ~f'(OdC  if‘(C)dC  j 

Equation  (Al)  is  discretized  by  central  differences  on  a  rectangular  grid  having  uniform 
^pacing  in  (r.C),  corresponding  to  variable  spacing  in  z.  The  finite-dimensional  version  of 
equation  (Al )  is — 


where: 


DO  +  il(0,A)  =  0 


<<£);,  -  <t>(  r,.  ) 


'ill..  H  \  ili'u  -  '0>l(  =  A  Pf(£Dyn i  -  r,-  Q((0),i.) 


D  =  l)r  +  \Y 


-(Tt:) 

ird  '-*) : 


j=q  and  i=p-  i 


j=q  and  i=p 


j=q  and  i=p+ 1 


otherwise 


!=p  j=q-i 
i=p  j=q 

i=P  J=9+1 
otherwise 

Dr.  DA  are  directional  operators  containing  entnes  relevant  to  the  denvanves  in  r.  _ 
di  jwtions.  Each  contains  three  nonzero  diagonals  in  a  block  structure  as  shown  in  figure 
A 1.  The  separation  or  D  into  us  directional  components  and  the  structure  of  these 
components  wiil  be  used  in  section  A. 2  below  to  form  efficient  and  consistent  concepts  of 
separation  of  variables  and  inner  product  for  the  discrete  problem  A2. 

<  FIGURE  A 1  ABOUT  HERE) 

The  contribution  of  boundary  points  with  Dirichlet  type  boundary  condiuon  is  placed 
in  12  when  non-zero:  for  Neumann  type  ooundarv  condition,  an  external  node  is  defined 

dO 

outside  the  boundarv.  its  value  uiven  bv  the  tirst  derivative  at  the  boundarv:  tor  — : — =»>.  we 
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set  On„i=Oi,.i  in  the  equation  for  the  boundary  node  .  The  singularity  of  (Al)  at  r=0  is 

not  explicitly  present  in  the  numerical  problem,  since  d>=0  at  r=0.  and  the  discretizauon  of 

( A I )  takes  place  only  in  the  interior  of  the  domain — a  finite  distance  from  r=0  (more  on  this 

singularity  in  section  A. 4). 

A. 2  Algebraic  Treatment  of  the  Bifurcation  and  Continuation 

In  the  neighborhood  of  any  bifurcation  point  A=(i.  we  expand  the  solution  to  i  A2)  in 
a  smaii  parameter  e.  as  in  <  1 8): 

<b  =  r  Q0  *  e-  Qi  -  . .  A  f  i 

A  =  U.  +■  £  Ki i  -  t: ~  K;  f  ... 

Substituting  into  ( A2)  and  coilecung  terms  of  like  power  in  e.  we  obtain  the 
sequence  of  equations,  anaiogous  to  ( 20).  the  first  three  of  which  are: 
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!  D  + 

u  E  :i  -U  ' 

'  !Qo=  !  D  +  s,n  |  Oo  =  n 

■  A  4a  i 

L 

(^1  =  -  K'n  P 

1  -  S' Oo  <£to 

'  A  4  b ) 

L  fi:  =  - 

■<iEvl 

*0,)  -  K'nP  ‘  A 

l  -  -S'-’doSi  -  S  ,(2oSoi2()  - 

A  Ac  i 

NIlP  ~ 

’Oofin 

where  the  derivatives  of  12  are  defined  following  1 19): 
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(0=0)  ,  ... 


Stk)  =  j x  P'k)  -  Q'k) 

The  eigenvalue  problem  (A4ai  can  be  solved  by  separation  of  variables:  let  Qo  have 
the  form,  euuivaient  to  (21) 


<£0)ij  —  '2Lr))i 


r  z 

’.vnere  %  and  may  be  called  "directional''  components  of  Qo  •  Substituting  into  t  A4a)  ana 
using  the  repeated  block  structure  of  D  we  obtain: 

tt(— %ill  "r  M-<E  ^plU  HSidl  J  <Q/))q 

(^))p{(DOlqlj(^)J)=0 

which  is  an  equality  of  two  rank- 1  matrices:  there  exists  therefore  a  scalar  3  such  that: 

*3  fO<))q  _r  *.QZ)|qij  (2n)j  =  ^  (A  5a) 


.ana: 


Ox-1 


.  .Dl 


'Dr)nl]1  -  ti  Lit.  -  .cr'Vi,,,  nidi  -  u  (Qi,)p  =  o 


A  5b) 


Note  that  P  1  is  invernble  since  it  is  a  diagonal  matrix  having  p(rj)  entries  on  its  diagonal, 
w  nich  are  positive. 

The  eigenvalues  3n>  M-nm  and  their  eigenvectors  are  found  from  <A5  )  by  the 
:ndiagonal  set  of  subroutines  from  Eispack. 

Detine  an  inner  product  as  follows: 

<x,x>  =  !rSy  A  ()i 

a  nere: 


'  SJ py  1 1  —  lS  )p,  (_S/  )lll 


and  Sr,  S/  are  diagonal:  Sr  =  tiiauj  .  Sz  =  ;iiag[f’(T,;]  (the  first  ana  last  eiements  may 

be  different  due  to  the  boundary  conditions);  both  are  positive  definite.  Under  this  inner 
product.  (Pa  Vl  is  seif-adioint  and  the  eigenvectors  of  <A4a)  are  onhogonai. 

Multiply  (A4b)  by  (P'1,VI  to  obtain  the  seif-adioint  form  of  the  operator,  and  we  can 
impose  a  solvability  condition: 

<  (  P'1^)  1  IS'*-'  6n  <JVq  +  ten  P' 1 J  OqJ  ,  >  =  0  (A  7) 


This  is  used  to  find  Kn.  as  in  (23).  To  soive  ( A4b)  for  £i ,  we  use  an  eigenvector 
expansion  method.  Similarly.  K)  and  Qican  be  found  from  (A4c).  and  so  on  to  any  desired 
order.  We  used  two  or  three  terms  of  (A3),  depending  on  the  type  ot  bifurcation 
encountered. 


r)d> 

Let  d>  = - .  and  take  the  derivadve  or  ( A2)  with  respect  to  A: 

r)  A 


LO  =  -PfO)  i  AS) 

This  is  a  system  of  ODE's  for  <|>(A),  equivalent  to  (29b):  given  initial  conditions 
lb(Apj  from  (A3)  or  the  weakiv-noniinear  esamate  (30),  this  linear  system  can  be  solved 
for  O.  We  use  Sparspak  (George  et  al.  1980)  to  solve  (A8)  and  a  Runge-Kutta  integrator 
(Press  et  al.  1986)  to  integrate  over  A. 

Given  d>0  and  A,  an  approximate  solution  of  (A2),  we  may  define  a  Newton  s 
Method  iteration  (Dennis  &  Schnabel.  1983]: 

=  S&ic  -  [D+AE^Q^r1  [D  Ok+GtOk  , A)  1  (A  9) 


w'here: 

a(£fi)ij 

■of’’™  =  '^'4=^ 

Under  some  mild  assumptions,  in  particular  that  <t>o  be  close  enough  to  the  exact 
■solution  and  that  the  Jacobian  oe  nonstngular  and  Dipschitz  connnuous.  this  iteranon  wiil 
converge  to  the  exact  soluuon  at  quadratic  rate  |  Dennis  &.  Schnabel,  theorem  5.2. 1  ].  If  we 
make  the  initial  estimate  (A3)  and  the  integraaon  (A8)  accurate  enough  and  stay  away  from 
bifurcation  points  (where  the  Jacobian  is  singular),  then  convergence  of  this  step  is 
practically  guaranteed. 

We  use  Sparspak  to  soive  each  step  of  ( A9).  and  a  line  search  algorithm  to  improve 
global  convergence  propernes.  The  structure  of  the  Jacobian  in  ( A9)  is  the  same  as  that  or  U 
in  (AX),  so  the  most  time-consuming  part  of  the  Sparspak  algorithm — -the  structure 
decomposition — needs  to  be  done  oniv  once. 


•n- 


A. 3  Numerical  Errors  and  Conversance 

Two  types  of  numerical  errors  need  to  be  considered:  the  discretization  error  (the 
difference  between  the  exact  analytic  solution  and  the  exact  solution  of  the  discretized 
system),  and  the  convergence  error  (the  difference  between  the  numbers  actually  obtained 
and  the  exact  solution  of  the  discretized  problem).  Convergence  errors  are  important  in  the 
corrector  step  only  since  those  occurring  in  the  predictor  steps — initial  estimate  and 
integration — are  irrelevant  when  the  corrector  step  converges. 

The  initial  esumate  need  not  be  very  accurate,  as  explained  above.  However,  if  it  is 
too  inaccurate,  then  the  corrector  procedure  may  converge  to  a  different  branch  or  not 
converge  to  a  solution  at  all.  The  initial  estimate  will  improved  as  £— D  in  (A3),  but  the 
integration  and  corrector  steps  wiil  lose  accuracy  as  the  Jacobian  becomes  singular  near  the 
bifurcation  point.  We  found  that  with  £=0.01.  taking  up  to  3  terms  of  the  senes  (A3)  leads 
to  convergence  of  ( A9)  and  reasonable  accuracy  for  the  subsequent  integration. 

The  convergence  errors  are  determined  by  the  stopping  criterion  of  the  corrector  step 

i  Ay): 


II  D  O  +  Q(O.A)  Ilf  <  52 


The  change  in  ilOII  in  the  last  Newton  step  is  usually  considered  to  be  of  the  same  order  as 
the  convergence  error.  For  5=1  O'8,  the  typical  value  was:  II  AO  11=1 0~5.  As  shown  below, 
this  is  much  smaller  than  the  estimate  for  the  discretization  error,  and  we  may  therefore  treat 
the  computed  solutions  as  "exact’’  solutions  of  the  discretized  problem. 

Equation  ( A2)  was  derived  using  central  differences,  and  is  second-order  correct  in 
Ar.  AC.  The  grid  transformations  z— >C  used  have  a  finite  derivative  everywhere,  and 
therefore  do  not  effect  the  order  of  the  truncation  error  (Kalnav  de  Rivas.  1972);  therefore  it 
is  considered  second-order  also  in  Az.  Higher-order  accuracy,  as  well  as  an  estimate  for  the 
discretization  error,  may  be  obtained  bi^ichardsoris  Extrapolation.  Computation  was 
repeated  for  sample  cases  with  grids  having  resolutions  of  (N^N,),  where  N.  &.  N,  e 
1 10.20,40):  9  different  grids,  followed  by  a  2-dimensional  extrapolation  to  —  0. 

Typical  values  of  the  relative  difference  of  the  second-order  solution  from  the  extrapolated 
results  are  presented  in  figure  A2:  these  differences  serve  as  an  estimate  for  the 
discretization  error. 


( FIGURE  A 2  ABOUT  HERE) 

To  further  vaiidate  the  above  estimate  of  the  discretization  error,  we  applied  the 
numencal  algorithm  to  a  problem  having  a  known  solution.  The  nonlinear  function  L2  in 
equation  ( i 6)  is  chosen  to  be: 

<t>2 


f2(0.r:A)  =  AO  + 


r  J  i  ( a  r  ) 


'.vnere  J-,  is  a  Bessei  function  or' order  i  ana  a  is  its  first  zero  i a  =5.8317...,.  Equation 
:  16)  with  this  LI  has  an  analytic  solution  which  is  qualitatively  similar  to  the  computed  tana 
weakly  nonlinear)  solitary  wuves: 

Oir.z;  =T-(£i  2-A)  r  J](a  ri  sech"(^T  z  V  a  :-Aj 

This  soiurion  bifurcates  from  the  trivial  <t>=0  branch  at  A=a  2  and  increases  in  magnitude  as 
A— >0.  (This  problem  does  not  necessaniy  correspond  to  a  physical  primary  tlow.  i  The 
discretization  errors  for  O.  Pot  and  the  first  two  axial  wavenumbers  are  presented  in 
figure  A3. 

( FIGURE  A3  ABOUT  HERE) 


The  error  in  the  perturbation  streamfunction  with  the  20x20  grid  is  close  to  2%  for 
the  test  function  and  less  than  IT  for  the  Richardson-extrapolated  case.  For  velocities 
(computed  from  the  streamfunction  by  central  differences),  the  discretization  error  is  larger, 
but  still  not  exceeding  a  few  percent  on  a  20x20  grid.  The  results  in  figure  A2(b)  are  for  the 
axial  velocity  w=du//r3r,  which  is  strongly  dependent  on  r-resolution.  The  errors  in  the 

values  of  the  axial  wavenumbers  V&i  and  the  bifurcation  points  p.^  are  similarly  of  order 
1  %  for  the  same  level  of  resolution.  The  20x20  grid  was  therefore  the  standard  in  most  of 
our  computations. 


A.4  The  Singularity  at  r=0  -H? 

-  r  J 

Equation  ( 16)  has  a  singularity  on  the  axis  r=0.  and  construction  of  a  numerical 
scheme  as  well  as  interpretation  of  the  results  should  take  that  singularity  into  account.  The 
discretization  ( A2)  makes  explicit  use  of  the  boundary  condition  (16c)  at  r=0.  and  applies 
'AD  only  to  intenor  gnd  points:  the  singularity  is  thus  avoided.  However,  if  the  end  is 

refined  unnl  j  becomes  very  iarge  at  the  first  gnd  point  off  the  axis,  then  the  matnccs 

involved  will  become  unbalanced  and  numencal  accuracy  will  detenorate.  1".  our  case  such 
fine  gnds  are  not  necessary'  since  Richarason  s  extrapolation  seems  to  show  convergence 
before  very  large  numbers  occur. 

The  singularity  is  encountered  again  when  we  compute  the  axial  velocity  at  the  axis, 
which  is  used  as  a  measure  of  the  perturbation  size.  The  definition:  w=r ^chiz/dr  cannot  be 
applied  directly  at  r=0:  two  numencal  schemes  are  used,  and  the  values  obtained  for  wt  0.0) 
agree  to  within  a  few  percent. 

The  first  method  is  a  quadratic  extrapolation  of  w  values  from  interior  gnd  points, 
coupled  with  the  condition:  cnv/dr=0  at  r=<).  The  quadratic  function  satisfying  tiiis  condition 
and  passing  through  the  first  two  interior  gnd  points  is: 


wm.Av,^n[4-(^):]+yv,24r)[(C-)--!] 


leadinti  to: 


-wt  Ar)  -  '.v(2Ar) 

wtO)  = - - - 

■> 

The  second  metnod  applies  Stokes'  theorem  to  a  rectangular  loop  of  dimensions  tor. 
6z)  touching  the  axis  r=0  and  centered  about  the  line  of  symmetry  z=0.  To  reduce  me  error 
associated  with  numerical  integration  over  a  finite  rectangle,  we  let  6 z—a)  and  obtain  a 
balance  involving  r-integration  only.  Tne  vomcity  integral  can  be  expanded  in  powers  of 
oz.  T)  being  the  azimuthal  component  of  the  voracity: 

5r  <34/2 

|oo»ds  =  |  |  j"  r)(r,0)+z:r^r,0)+...]  dzdr 

s'  0  -0/72 L  <!L  J 

5r 

=  5z  i  r\(r,0)  dr  +  Olbz^i 

o' 

A  similar  expansion  is  done  for  the  circulation  integrals: 


oz/2 

j3  u»dl  =  J[w(0.z)-w(5r.z)]  dz 

-oz/2 

or 

+  J"[u(r,bz/2)-u(r.-Sz/2)l  dr 
0 


=  5z 


w( 0,0)  -  wtdr.O)  - 


5r 

ri&<f 


rdz 


tU,0)  dr 


+  OfSz^i 


Companng  the  leading  terms  in  oz.  we  obtain: 


or 


d2\l/ 


wtO.O)  =  wtdr.O)  -t-  j,  T](r,0)  +  — -^nr.O) 


r  rJz^ 


dr 


The  integral  was  computed  using  the  Simpson  T-rule.  and  the  two  expressions  for 

wtO.O)  were  compared  for  or=Ar  and  2Ar.  The  differences  were  of  order  1  'T  in  most 
cases,  and  increased  up  to  5nc  oniv  as  A-^u  (where  the  perturbation  is  small  and  roundoff 
error  becomes  significant)  and  as  A— ?()  (where  large  radial  gradients  require  increased 
resolution).  We  therefore  used  the  sunnier  quadratic  extrapolation  form  mroughout.  This 
comparison  aiso  serves  as  an  additional  check  on  the  convergence  of  the  numerical  results 
near  tne  singular  line  r=0. 


:4 
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igure  Captions 


Figure  i.  Bifurcation  diagram  tor  the  first  two  columnar  Drancnes.The  extreme 
perturbation  axial  velocity  occurs  away  from  the  axis  on  part  of  the  second 
branch. 

Figure  2.  Velocities  on  the  second  columnar  branch,  for  different  vaiues  of  the  swiri 
parameter  A.  (a)  axial  velocity'  tb)  azimuthal  velocity. 

Figure  3.  Velocities  on  the  supercritical  side  of  the  principal  columnar  branch,  (a) 
axial  velocity  (b)  azimuthal  velocity. 

Figure  4.  Velocities  on  the  subcritical  side  of  the  principal  columnar  branch,  (ai  axial 
velocity  tb)  azimuthal  velocity. 

Figure  5.  Exchange  of  criticality  between  the  primary  flow  i  — o — >  and  the  principal 
conjugate  branch  ( — * — -.  using  Benjamin's  Froude  number  criterion. 

Figure  6.  Comparison  of  computed  (  o  )  and  asvmptonc  t  —  >  solutions  in  the  limit 

A— >0  on  the  principal  branch,  (a)  location  of  the  boundary  layer  tb)  maximum 
of  the  streamfunction  u i*  tc)  layer  thickness  6. 

Figure  7.  Comparison  of  computed  and  asymptotic  velocity  profiles  for  the  principal 
branch,  at  A=().()01|loo  . 

Figure  8.  Bifurcation  diagram  for  the  periodic  branches  of  Uni  and  u(>2  :  u=2  for 

better  separation  of  the  bifurcation  points. - principal  columnar  branch  : 

— o —  L=  10  wave  centered  at  z=L/2:  — •—  L=  10  wave  centered  at  z=U: 

— Q —  L=5  wave  centered  at  z=L/2:  — ■ —  L=5  wave  centered  at  z=(). 

Figure  9.  Approach  of  the  periodic  solutions  to  localized  w'aves.  <ai  fixed  lengtn  L=o: 
wave  trough  becomes  localized  as  the  wave  amplitude  increases,  tb)  fixed 
amplitude  A/)i=1.10:  the  half-height-length  of  the  wave  approaches  a  constant 
value  independent  of  the  computational  domain  length. 

Figure  10.  Meridional  streamlines  of  periodic  solutions:  L=3.  a=14.  contour 

intervals  ot  0.1  .  iai  (A-ui/(j  =  1.5  .  just  before  the  appearance  of  stagnation 
points,  lb)  (A— (J.)/p  =  l  ,b.  a  recirculation  bubble  appears,  to  (A-ui/(i=2.5.  the 
recirculation  region  grow  s  with  the  perturbation  amplitude,  tdi  detail  of  the 
bubble  in  to. 


Figure  il.  Velocities  on  me  periodic  nrancn  of  u,r,  ai  me  wave  center  z=u.  <ai  axiai 

velocity  tb)  azimuthal  velocity  to  dii’ference  in  axial  velocity  between  the  wave 
center  and  tab. 

Figure  12.  Solitary'  wave  behaviour  oh  the  periodic  solutions:  the  dependence  or 

effective  wave  length  t half-height  length)  on  the  wave  amplitude,  compared  to 
the -0.5  slope  for  an  exact  solitary' wave,  o  L=b:  •  L=20:  +  L=l()().  (a) 
wave  amplitude  measured  by  extreme  axial  velocity  tb)  wave  amplitude 
measured  by  the  distance  from  the  bifurcation  point. 

Figure  13.  Convergence  of  the  wave  crests  of  the  periodic  solution  to  the  columnar 
solution  of  the  same  amplitude,  t  a )  fixed  length  L=6,  the  difference  vanishes  as 
the  wave  amplitude  increases,  (b)  Fixed  amplitude  (A-poi Vdoi-0.05.  as  the 
computational  domain  iength  increases. 

Figure  14.  Bifurcation  diagram,  solitary'  wave  branch.  The  extreme  perturbation  axial 
velocity  is  off  the  axis  when  A<0.16uno  . 

Figure  15.  Axial  velocity  at  wave  center  z=0  on  the  solitary  wave  branch 

Figure  16.  Meridional  streamlines  of  solitary  wave  solutions.  L=5.  contour  intervals 
0.05.  (a;  A/fioo  =0 .80  (b;  A/doo  =0.70.  a  small  recirculation  bubble  appears 
(c)  A/p(x)  =0.001,  a  large  recirculation  bubble,  (d)  detail  of  the  bubble  in  ib). 
(e)  detail  of  the  bubble  in  <c). 

Figure  17.  Comparison  of  solitary'  wave  amplitude  of  the  computed  (  o  )  and  weakly- 
nonlinear  t — 1  solutions,  i ai  bifurcation  diagram  ( hi  detail  showing  the 
bifurcation  and  the  amplitudes  of  flow  reversal,  (o  half-height-length  vs. 
amplitude  ot  the  computed  wave,  compared  to  the  -0.5  slope  of  the  exact 
solitary'  wave. 

Figure  18.  Comparison  of  perturbation  streamlines  of  the  computed  ( —  t  vs.  the 

weakly-nonlinear  < . i  solitary  waves,  (a)  A/doo=0.90.  contour  intervals 

0.004  <b)  A/u<x)=0. 50.  contour  intervals  0.012  (c>  A/doo=0.20.  contour 
intervals  0.02  . 

Figure  AI.  Structure  of  the  matrices  tat  D‘  ,  <b)  D‘ 

:XllJ  ■■  M’J 


Figure  A2.  Variation  with  mesn  size  of  errors  relative  to  Richardson-extrapolated 

values  at  A/Uoo=O.X0:  •  actual  computation  with  this  mesh:  o  extrapolated, 
i a)  error  in  0(0.3. 0)  (hi  error  in  wit). 0)  io  error  in  Uno  . 


ad- 


Figure  A 3.  Variation  with  mesh  size  or  errors  relative  to  tne  exact  solution  or  the  test 
problem:  •  actual  computation  with  this  mesh:  o  extrapolated,  tai  maximum 
error  in  <t>(r.zi  at  A Ja  :=0.S0  (hi  error  in  axiai  wavenumbers  icj  error  in  Uui 

(d)  error  in  Uoi  • 
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A  method  for  computing  leading  eigenvalues  (having  the  largest  real  pan)  and  their 
eigenvectors  for  large  generalized  eigenvalue  problems  is  presented.  A  linear  fractional 
transformation  is  used  to  map  a  group  of  leading  eigenvalues  into  dominant  eigenvalues 
(having  the  largest  modulus).  The  Dominant  eigenvalues  of  the  transformed  problem  are 
computed  by  Stewan’s  (1976)  Simultaneous  Iteration.  Each  iteration  involves  matrix- 
vector  multiplication  and  the  solution  of  a  linear  system,  which  can  be  done  efficiently  if  the 
matrices  involved  are  sparse  or  have  some  special  structure.  Convergence  properties  are 
similar  to  those  of  the  inverse  power  iteration:  the  method  requires  an  estimate  for  the 
region  in  the  complex  plane  containing  the  desired  eigenvalues,  and  converges  rapidly 
when  a  good  estimate  is  available.  The  amount  of  work  is  also  comparable  to  that  of  the 
basic  inverse  iteration,  which  is  significantly  less  than  that  required  for  full  eigensolution. 
Examples  from  hydrodynamic  stability  demonstrate  convergence  rates,  computation  time 
and  the  ability  to  resolve  simultaneously  groups  of  leading  eigenvalues. 


Computation  of  Leading  Eigenspaces  for 
Generalized  Eigenvalue  Problems 


1.  Introduction 

A  generalized  eigenvalue  problem  has  the  form: 

A  x  =  c  B  x  (l.l) 

where  A,  B  €  Gnxn  are  general  complex  matrices.  In  many  applications  these  matrices  will 
have  some  useful  structure,  such  as  symmetry  or  sparsity. 

Let  the  Leading  Eigenvalues  of  (1.1)  be  those  having  the  largest  real  part  ;  the  more 
common  term.  Dominant  Eigenvalues,  refers  to  those  having  the  largest  modulus  .  In 
some  applications,  only  a  few  leading  eigenvalues  of  (1.1)  are  sought;  for  example,  in 
hydrodynamic  stability  problems,  the  real  pan  of  o  is  the  growth  rate,  and  the  eigenvectors 
of  the  leading  eigenvalues  represent  the  most  unstable  modes. 

Traditional  methods  for  solving  (1.1)  usually  involve  finding  all  the  eigenvalues,  using 
the  Q-Z  algorithm  (see  IMSL  or  other  numerical  analysis  libraries)  and  then  soning  by  the 
real  pan.  This  involves  0(n3)  work,  where  n  is  the  order  of  the  matrices,  and  becomes 
expensive  or  impractical  for  large  n;  little  or  no  advantage  can  be  taken  of  sparsity  or  other 
structure  of  A  and  B. 

Several  methods  exist  for  extracting  selected  eigenvalues  and  eigenvectors  of  standard 
eigenvalue  problems,  i.e.  when  B  is  invertible  (see,  for  example,  Golub  and  Van  Loan, 
1983.)  Power  and  Lanczos  methods  compute  the  dominant  eigenvalues;  inverse  iteration 
can  find  the  eigenvalues  closest  to  a  given  point  in  the  complex  plane  and  their 
eigenvectors.  These  are  not  directly  applicable  to  the  problem  of  computing  the  leading 
eigenvalues. 

Recently,  an  integration  method  was  proposed  (Goldhirsch  et  al.  1987)  for  the  leading 
eigenvalues  of  a  standard  eigenvalue  problem  (where  B  is  invertible.)  This  method  is 
simple  and  elegant;  however,  its  convergence  may  become  very  slow  (or,  alternatively,  the 
size  of  the  reduced  problem  may  become  too  large)  if  the  separation  of  the  eigenvalues  is 
small.  Another  problem  may  arise  if  the  problem  is  defective,  i.e.  a  leading  eigenvalue  has 
generalized  eigenvectors;  in  this  case,  the  integration  method  may  return  inconclusive  or 
inconsistent  results. 


2.  The  Dominance  Mapping  Method 

This  method  attempts  to  address  the  problems  (1.1)  which  are  not  solved  efficiendy  by 
the  other  methods  mentioned  above.  It  will  work  for  singular  A  and  B;  for  defective 
problems;  it  will  take  full  advantage  of  the  structure  of  the  matrices;  and  it  allows  some 
control  over  convergence  rates.  There  are  a  few  restrictions,  however,  which  will  be 
discussed  below. 

The  eigenvalues  in  the  complex  a-plane  can  be  mapped  to  a  H-plane  by  the  linear 
fractional  transformation: 


o  =  B  +  a - 

X+l 


(2.1) 


where  a  is  a  real  positive,  and  p  a  complex,  constant.  The  important  effect  of  this  linear 
fractional  mapping  is  to  map  the  half-plane  to  the  left  of  a=p  to  the  inside  of  the  unit  circle 
in  the  X  plane,  as  seen  in  fig.  1.  If  m  leading  eigenvalues  are  required,  arid  we  select  P 
such  that: 


Re(<li)  (  >  ^e<^) 

{  <  Re(P)  i -m  +1  ...n 

then  the  corresponding  m  eigenvalues  will  be  dominant  in  the  X  plane: 


Figure  1:  The  Dominance  Mapping  (2.1) 


The  eigenvalue  problem  for  X  is  in  standard  form: 

C  u  =  Cf 1  C2  u  =  X  u  (2.2) 
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where: 


Ci  =  -  [A-«x  +(3  )  B] 

C2  =  [A+(a  -13 )  B] 

The  problem  of  computing  the  leading  eigenvalues  of  (1.1)  becomes  that  of  computing 
the  dominant  eigenvalues  of  (2.2);  the  methods  mentioned  in  §1  can  now  be  applied.  We 
used  Stewart’s  (1976)  version  of  Simultaneous  Iteration,  which  applies  to  the  most 
general,  non-hermitian  C. 

A  transformation  similar  to  (2.1)  was  proposed  by  Jennings  (1977),  in  the  context  of 
convening  a  quadratic  eigenvalue  problem  to  standard  form.  Jennings  (and  no  one  else,  to 
the  best  of  the  author’s  knowledge)  has  not  made  the  second  step  of  applying  a  dominant 
eigenvalue  method  to  a  transformed  problem  equivalent  to  (2.2). 

The  mapping  constants  a  and  [3  allow  the  user  some  control  over  the  rate  of 
convergence  and  the  order  in  which  the  leading  eigenvalues  emerge  during  the  iteration. 
The  user  must  have  an  estimate  of  where  in  the  complex  plane  the  leading  eigenvalues 
reside;  [3  is  set  to  the  left  of  this  region.  The  point  c  =  [3+a  is  a  singular  point  of  (2.1) 
which  maps  to  infinity  in  the  X  plane;  eigenvalues  close  to  c  will  map  to  very  large 
modulus  in  the  X  plane,  and  will  converge  rapidly  during  the  iteration  of  (2.2).  a  should  be 
set,  therefore,  so  that  c  is  near  the  center  of  the  leading  region  or  near  the  most  important 
eigenvalue. 

The  following  algorithm  computes  m  leading  eigenvalues  and  eigenvectors  of  (1.1), 
using  the  Dominance  Mapping: 

1.  estimate  leading  region;  select  a,  [3 

2.  perform  L-U  decomposition  of  Ci  =  -  [A-(a  +|3)  B]; 

(use  the  structure  of  A  and  B  ) 

3.  select  m  initial  column  vectors  6  (Enxm 

4.  Simultaneous  Iteration  on  Cu  =  ilu: 
for  each  multiplication  u^k+^  =  C  do: 

4.1  multiply:  v  =  C2u^ 

4.2  solve  the  system:  u'"lc+^  =  v 

5.  map  converged  X\  — »  Gj . 

3.  Singularities  in  the  Dominance  Mapping 

The  algorithm  of  §2  may  fail  in  two  cases,  corresponding  to  the  two  singularities  of  the 
mapping  (2.1):  the  point  o=c  ,  which  maps  to  infinity  in  X,  and  X=-l  which  maps  to 
infinity  in  a. 
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When  I c  -a;  i  <  £  for  some  i  <m,  for  a  small  (machine-dependent)  £,  then  the  matrix  C 
will  be  ill-conditioned  or  numerically  singular.  This  is  easily  remedied  by  a  small  change  in 
a,  which  does  not  significantly  affect  any  other  properties  of  the  mapping. 

When  IlmfcJi  -  c  )  I  »  1  for  some  i  <m ,  the  corresponding  ^-eigenvalue  is  close  to  the 
singular  point  A.=-l.  This  implies  that  its  separation  from  the  subdominant  eigenvalues 
inside  the  unit  circle  is  small,  often  so  small  that  convergence  is  impractical.  Some 
improvement  may  result  if  we  increase  a;  but  this  may  decrease  the  modulus  of  other 
dominant  X-eigenvalues  and  slow  down  their  convergence.  In  a  case  where  leading 
eigenvalues  are  widely  separated  in  the  imaginary  direction,  it  may  be  necessary  therefore 
to  restart  the  iteration  with  different  (3  values  to  resoive  separate  clusters  of  leading 
eigenvalues. 

4.  Example: 

The  performance  of  the  DM  method  can  be  demonstrated  by  observing  the  amount  of 
work  needed  to  resolve  a  fixed  subset  of  leading  eigenvalues,  as  the  order  of  the  problem 
increases.  The  following  example  includes  tridiagonal  matrices  of  increasing  size  N,  all 
having  two  leading  eigenvalues: 

(7!  =  1.2  (4.1) 

<32  =  1-1 

Re(<Jj)  <1  for  i  =  2 ...n  . 

Selecting  a  =  0.3,  (3  =  1.0  isolates  aj,  <32  ■  The  problem  was  solved  first  using  the 
traditional  QZ  routines  (IMSL),  then  using  the  DM  method  but  treating  the  matrices  as 
dense,  and  finally  taking  full  advantage  of  the  structure.  The  results  are  shown  in  Figure  2. 

The  savings  in  computing  time  relative  to  the  full  eigensolution  can  be  significant:  at 
n  =100,  only  j  of  the  work  is  necessary  even  without  exploiting  the  band  structure;  the 

work  is  reduced  by  more  than  an  order  of  magnitude  when  the  structure  is  used. 
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Figure  2:  Comparison  of  work  to  resolve  2  leading  eigenvalues  of  (4.1) 


5.  Application  to  the  Orr-Sommerfeld  Equation 

The  Orr-Sommerfeld  equation: 

(D2-  a2)"  \\f  =  iaR  [(U-c)(D2-a2)\|/  -  U  y  ]  (5.1) 

describes  the  hydrodynamic  stability  of  parallel  shear  flow  (see,  for  example,  Drazin  and 
Reid  1981.)  High  accuracy  eigenvalues  were  computed  by  Orszag  (1971)  for  plane 
Poiseuille  flow  with  R=  10000  (Reynolds  number)  and  a=l  (streamwise  wavenumber.) 
The  location  of  the  twelve  least  stable  eigenvalues  are  shown  in  figure  3. 

Equation  (5.1)  is  discretized  using  central  differences  (a  spectral  method  may  be  more 
appropriate  in  this  specific  case,  as  in  Orszag  (1971),  but  the  banded  finite-difference 
matrix  is  a  good  example  of  candidate  problems  for  the  DM  method.)  The  eigenvalue  c  is 
replaced  by  G=-ic  ,  to  conform  with  the  definitions  of  (1.1). 

When  Im((3)=Im(Oi),  eigenvalues  1  and  4  were  the  first  to  converge;  2,  3,  5  and  6  took 
longer  ,  a  converge,  since  the  imaginary  part  separation  brought  their  A.  counterparts  close  to 
the  singular  point  X=-l.  For  Im((3)=Im(02),  the  order  was  reversed:  first  eigenvalues  2.  3. 
5  and  6  and  then  1  and  4.  In  both  cases,  the  first  group  converged  within  10  to  15 
iterations,  regardless  of  the  number  of  grid  points. 

The  error  associated  with  convergence  of  the  /.-iteration  was  not  significant  in  our 
computations.  Using  a  stopping  criterion  of  IlCu-Uull  <  10*  A  the  leading  a-eigenvalues 
were  converged  to  at  least  5  digits.  The  discretization  error  of  the  finite-differencing 
(compared  to  Orszag’s  results)  is  proportional  to  Ax^,  as  expected.  The  time  to  resolve  the 
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Figure  4:  (a)  the  time  to  compute  the  most  unstable  eigenvalue  of  (5.1 ) 
(b)  discretization  error  vs.  the  grid  interval  Ax 
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6.  Conclusion 

Using  the  Dominance  Mapping  and  a  Power  Iteration  method  we  can  compute  leading 
eigenvalues  and  eigenvectors  of  large  generalized  eigenvalue  problems.  This  method  can  be 
more  efficient  than  a  full  eigensolution  even  for  a  general  pair  of  matrices,  but  is  especially 
attractive  when  the  matrices  have  a  structure  that  can  save  work  in  the  Gaussian  elimination 
and  matrix  multiplication  steps.  The  DM  method  can  be  applied  to  singular  and  defective 
problems  that  may  cause  failure  or  slow  convergence  in  other  methods. 

Use  of  the  DM  method  is  restricted,  however,  to  cases  where  an  estimate  for  the 
leading  eigenvalues  is  available.  When  this  estimate  shows  a  wide  distribution  of  leading 
eigenvalues  along  the  imaginary  direction,  several  passes  may  be  necessary  with  different 
mapping  parameters  to  properly  resolve  all  leading  eigenvalues,  as  demonstrated  for  the 
Orr-Sommerfeld  problem. 
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